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Preface 


As  a  rcsalt  of  ioforaal  discassioas  darii^  tlie  1975  Methoiem  aad  Verfaluiea  der 
aathcnatiscliea  Physik  wetii^  at  Ibervolfadi  it  vas  a|ipareat  tkat  tkere  va$  aa 
outstanding  need  for  aa  iateraatioaal  scaiaar  gtoap  vkick  coaceracd  itself 
with  probleas  vkick  kad  been  attcapted  bat  vkick  still  stabboraly  resisted 
solutions  rather  thaa  vitk  probleas  vhidk  kad  beea  solved. 

The  first  suck  scaiaar  vas  keld  at  tke  Ifahersity  of  Bona  ia  Septeidier 
1976  under  the  heading  of  lived*  Boaadary  Valve  Probleas.  Since  then  a  nnaber 
of  siailar  neetings  on  a  variety  of  different  topics  have  been  held  at  either 
the  Ifniversity  of  Bonn  or  the  University  of  Strathclyde.  There  is  nov  every 
indication  that  siailar  aeetings  vill  be  held  r^larly  aad,  aoreover,  that 
they  will  bccoae  increasingly  anltidisciplinary.  Consequently,  it  has  been 
decided  to  aake  the  vork  of  these  aeetings  aore  readily  ai’ailable  by 
publishing  proceedings.  This  volnae  coaprises  the  proceedings  of  the  aost 
recent  aeeting  held  at  Boss  Priory,  Ihiiversity  of  Strathclyde.  Ve  are  aost 
grateful  to  Longnan  Group  United  for  agreeing  to  publish  these  proceedings. 

These  various  aeetings  could  not  take  place  without  a  considerable  anount 
of  support  and  assistance  fron  a  nnaber  of  sources.  In  this  connection  ve 
are  particularly  grateful  to 

The  Royal  Society  of  Edinburgh 
The  Edinburgh  Matheaatical  Society 
The  U-S.  Amy 


for  tlieir  fiaouscial  sapport.  Ikuiks:  are  also  die  to  all  tkose  cclleagies  vko 
helped  vitli  tiie  refereeiig  of  tke  p^ers.  last  bit  bv  so  acais  least,  vc 
record  oir  appreciatioi  of  tbe  vorfc  doie  Inr  tbe  sccxetarjal  staff  of  tbc 
Departaeit  of  Xatheaatics  ii  tbe  ITiiversity  of  Stratbclyde,  especially  Xrs. 
Mary  Sei^eait  vbose  qiiet  efficieicy  helped  ii  so  my  vays  to  ensorc  the 
siootb  niiiig  of  the  coifereice. 

tfiiversity  of  Strathclyde  G.  F.  Eoach 

Clasgok'  G1  IXH 
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W.R.  BRECKON 

Measurement  and  reconstruction  in 
electrical  impedance  tomography 

INTIODUCriON 

In  this  paper  we  consider  the  inverse  problea  of  Electrical  lapedance 
Tonography  (EXT).  In  this  aedical  inaging  technique  a  pattern  of  current  is 
applied  to  the  surface  dtt  of  the  body  fl  and  neasurenents  of -the  resulting 
electric  potential  are  made  on  A  survey  of  clinical  application  of  this 
technique  can  be  found  in  Brown,  Barber  and  Seagar  [1]  and  comprehensive 
infomation  on  various  aspects  of  the  technique  in  the  proceedings  [2]  [3] . 

Mathematically  the  problem  can  be  formulated  as  follows.  Let  7  e  L“(ft) 
be  the  electrical  conductivity,  which  satisfies  7(x)  >  c  >  0  for  (almost)  all 
X  €  Q.  The  potential  u  then  satisfies 

V.7VU  =  0 

being  a  combination  of  Ohm's  law  and  Kirchoff's  current  law.  The  current 
density  j  on  ^  is  given  by 

j  =  -tV„» 

where  n  is  the  outward  unit  normal  on  30.  Vhere  convenient  we  will  use  the 
operator  =  V.7V.  To  solve  the  equation  L^u  =  0  it  is  sufficient  either  to 
specify  the  Neumann  condition  j  (together  with  an  additional  condition  on  the 
potential,  such  as  u(p)  =  0  for  some  p)  or  the  Dirichlet  conditions  u|30. 

When  one  of  the  sufficient  conditions  is  specified  the  other,  or 
complimentary,  boundary  condition  is  determined  also  for  a  given  7. 

O  O 

Whilst  ft  is  clearly  a  domain  in  R  the  case  of  ftcR  is  also  considered. 

The  question  of  the  possibility  of  identification  of  7  from  a  knowledge  of  all 
pairs  (j,u|3ft),  that  is  the  uniqueness  of  solution  of  the  inverse  problems, 
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has  been  considered  by  Kohn  and  Vogelias  [4] ,  who  proved  uniqueness  for  a 
piecewise  analytic  conductivity  and  Sylvester  and  Uilnan  [5]  who  proved  the 
sane  for  yc 

Ve  will  consider  two  aspects  of  the  inverse  problen  here  Ve  present  a 
suggestion  for  a  systea  of  optiaal  aeasureaents  and  a  linearisation  approach 
to  the  solution  of  the  inverse  problem.  In  the  latter  we  present  an  example 
to  show  that  the  linearised  inverse  problea  is  not  the  inversion  of  a 
Generalised  Radon  Transform  as  some  authors  have  assumed. 

Two- norm  optiaal  measurement 

The  ideal  case  considered  by  Kohn  and  Vogelius  and  by  Sylvester  and  Uhlman 
assume  that  we  have  a  perfect  set  of  measurements,  that  is  all  pairs  (j,u|dft) 
are  known.  This  is  equivalent  to  knowing  the  transfer  impedance  operator 

where  W  we  consider  the  finite  energy  case  u  e  then 

is  a  pseudo  differential  operator  which  is  compact  and 
self  adjdoint  as  a  map  -*  The  mapping  7  is  non-linear  and 

it  is  this  mapping  we  seek  to  invert. 

In  practice  we  can  only  apply  a  finite  number  of  current  patterns  and 
measure  voltage  only  at  a  finite  number  of  points.  The  question  arises: 
which  are  the  best  measurements  to  make? 

Isaacson  [6]  gives  one  answer  to  this  question  which  he  frames  in  terms  of 
a  measure  he  calls  distinguishability.  If  7I  and  72  are  two  conductivities 
then  they  are  distinguishable  by  measurements  of  precision  e  if  there  is  a 
current  density  jG  for  which 

^(j)  =  l|R.yij-R.y2jll/Njll>e 

the  number  ^(j)  is  called  the  distinguishability.  The  best  currents  in  the 
sense  of  Isaacson  are  those  which  maximise  ^(j) 
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where 


^(j)  =  snpMR^jj-Il^jll 

Mil  1=1 

=  sap<j,D^j> 

I  111  1=1 

»  = 

The  aap  D  is  a  coMpact,  self  adjoint  pseudo  differential  operator 
H®(dQ)  -•  H®(^) .  It  has  a  conplete  set  of  orthonormal  eigen  functions 
>  with  eigen  values  -♦  0  as  k  ^  m.  From  the 

;  min- max  principle  one  can  deduce  that  the  largest  distinguishability  possible 

;  is  A^  which  is  achieved  when  j  is  an  eigen  function  with  this  eigen  value. 

Isaacson’s  algorithm  for  calculating  this  optimal  current  is  based  on  the 
power  method  (see  for  example  [7]).  We  will  take  7I  to  be  the  (unknown) 
conductivity  of  the  body  and  ^2  as  the  best  available  guess  for  the 
j  conductivity.  The  method  is  an  iterative  process  which  involves  repeated 
I  measurement  and  can  be  expressed  as  follows: 

'  Guess  (where  =  1) 

i  Repeat 

i  Apply  j^^^  and  measure  v^”^  = 

Compute  v~^"^  = 

Set  A^  =  |lv(")-v~(")|| 

Set  =  v(")-v~("))/A„ 

Until 

The  reasoning  behind  making  repeated  measurements  of  v  rather  than  taking  a 
I  basis  of  currents  and  working  with  the  resulting  matrix  for  is  that  the 

I  measurement  process  involves  error  and  so  is  not  truly  linear.  However  the 

power  method  only  produces  the  largest  eigen  function,  which  would  give  only 
1  one  measurement  with  which  to  estimate  7I.  It  would  be  desirable  to  have  a 

j 
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basis  of  currents  at  least  spanning  the  sane  space  as  the  eigen  functionns 
with  Aj^<G.  For  this  reason  it  is  desirable  to  have  a  procedure  which 
computes  all  the  eigen  functions  jk  with  Aj^>e.  The  following  procedure  has 
been  found  to  work  in  computer  simulations  with  pseudo  random  noise  in  the 
measurement  stage. 

Guess  orthonormal  basis  with  [  jl  =  0) 


Repeat 

Measure  =  R^^jk^"\  compute  ^or  all  k 

Compute  rj^j^  =  <jk,vk("^-v~j^^"^> 

Calculate  the  eigen  system  for  R  =  {r^j^] ,  RU  =  Ua 
Set  =  S  Uyjj 

Until  1|R-A||<6 

In  numerical  experiments  little  was  gained  after  two  iterations. 

Point  optimal  currents 

Isaacson's  criterion  for  best  currents  gives  the  best  currents  only  in  the 
sense  of  optimising  the  two  norms  of  the  voltage  data  measured  for  a  given 
current  pattern.  In  practical  systems  each  measurement  of  voltage  is  made 
separately.  It  is  of  interest  therefore  to  find  a  current  pattern  which 
optimises  the  voltage  difference  at  a  point  on  Xl  between  bodies  of 
conductivity  'fl  and  72. 

Let  pG  do,  be  the  point  at  which  we  make  the  measurment 
v(p)-v~(p)  =  R.yij)(p)-(R^2j)(P)-  ^ 

rii)  =  lv(p)-v~(p)| 
is  optimised  subject  to  llj||  =  1  and  =  0. 

We  can  express  the  current  pattern  in  terms  of  the  eigen  functions  as 
j  =  Sjjjjj*  The  optimisation  problem  is  then 
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maxiiise  Saj^-^j^jjj(p) 

2 

subject  to  Eaj^  =  1. 

Using  Lagrangian  procedure  we  obtain 

\ 

where 

c  =  l/(5)(Akjk(p))^). 

The  measurement  procedure  would  then  be  first  to  calculate  the  eigen  functions 
using  the  method  of  the  previous  section  then  to  make  a  measurement  at  p  apply 
the  current 

j  =  S  CAj^jk(P)jk 

to  give  the  voltage  v(p)-v~(p). 

2 

In  the  simple  case  ft  =  {xeR  :|x|<l}  with  72  =  1, 

(<r  lxl<p 

7l(x)  = 

4  |x|>p 

the  eigen  values  can  be  explicitly  calculated  as 

where  p  =  (l-<r)/(l+<r)  and  the  eigen  functions  are  cos  k9  and  sin  k^. 

A  plot  of  the  point  optimal  current  for  this  case  with  p  =  0  and  various 
values  of  p  is  shown  in  Fig  1.  In  the  case  of  small  objects  in  the  centre 
this  procedure  deviates  least  from  that  of  Isaacson. 


Point  Optimal  Gurent  sigma  =  0.300  r  =  0.1  to  0.9. 

Figure  1 


Linearisation 

The  mapping  7  is  non-linear.  To  solve  the  inverse  problem  it  is 
desirable  to  linearise  this  mapping.  The  original  statement  of  the 
linearisation  appears  in  Calderon  [8].  In  this  chapter  Calderon's  techniques 
are  used  and  elaborated  upon  to  give  linearised  forms  of  the  forward  problem 
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both  in  its  direct  and  integral  form.  Since  the  original  reference  is  rather 
hard  to  find  it  is  hoped  that  this  will  help  readers  who  have  not  yet  located 
a  copy  of  this  legendary  Brazilian  conference  paper.  Calderon's  result  is 
extended  to  give  a  'Neumann  conditions  constant'  formulation.  All  published 
reconstruction  algorithms  rely  on  some  form  of  linearisation  and  yet  the 
approximations  used  are  not  always  justified.  With  this  in  mind  the  subject 
is  treated  in  some  detail  here.  We  find,  reassuringly,  that  they  are  all 
simply  the  Frechet  derivative  of  appropriately  defined  forward  mappings. 

If  the  conductivity  7  is  perturbed  to  7+^7  and  yet  one  form  of  sufficient 
boundary  data  for  u  is  kept  constant,  the  complimentary  boundary  data  will 
change.  For  example,  if  a  current  density  j  is  applied  resulting  in  a 
potential  u,  that  is  L^u  =  0,  with  -7VjjU  =  j  on  dfi,  then  when  7  is  changed  to 
7+^7,  the  potential  u  will  change  to  u+5u,  hence  L^^^^(u+5u)  =  0  with 
- (7+57)Vjj(u+^u)  =  j.  The  voltage  difference  on  the  boundary  <5uj^  will  be 
the  data  we  measure  in  an  attempt  to  detect  this  conductivity  change  so  we 
want  a  formula  for  5u  in  terms  of  ^7  (the  reverse  would  be  too  optimistic!), 
neglecting  higher  order  terms  in  ^7.  This  is  achieved  by  writing  5u  as  a 
series  in  Sj  and  truncating  after  the  linear  term.  This  series  involves  the 
linear  operators  which  depends  on  Sj  in  a  linear  way  and  which  is  the 
inverse  of  L^,(the  G  stands  for  Green's  function  of  course).  The 
non- specialist  reader,  if  daunted  by  manipulations  of  differential  and 
integral  operators  as  though  they  were  numbers,  may  like  to  think  of  them  as 
matrices,  as  they  would  be  if  we  passed  to  some  discrete  approximation  to  the 
operators. 

CHOICE  OF  SPACE  FOR  7 

Standard  elliptic  theory,  such  as  that  presented  in  Gilbarg  and  Trudinger 
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[[9],  requires  the  coefficients  of  the  partial  differential  equations  to  be  in 
L®(0)  which  is  a  fairly  weak  assumption.  To  guarantee  ellipticity  of  the 
operator  we  certainly  need  7  to  be  bounded  away  from  zero,  7  >  c  >  0 
(almost  everywhere).  In  the  following  results  we  will  need  to  be  able  to 
evaluate  |l7|dfi||,  that  is  we  need  to  estimate  the  magnitude  of  the 
conductivity  on  the  boundary.  In  L“(ft)  there  is  no  natural  restriction 
mapping  as  5Q  is  a  set  of  measure  zero.  While  L“(ft)  contains  extremely  nasty 
functions  it  has  an  extremely  strong  convergence  criterion.  We  would 
certainly  be  unwise  to  compare  images  on  the  basis  of  their  L“  distance. 

Natterer  [10]  suggests  that  the  appropriate  norms  with  which  to  compare 
two  dimensional  images  is  This  space  just  fails  to  include  the 

characteristic  functions  of  domains  with  sufficiently  regular  boundaries. 

The  weighting  of  high  frequency  terms  (or  if  you  like  the  inclusion  of  the 
l/2th  derivative)  weights  edges  more  strongly  than  the  simple  L  norm  and  this 
is  consistent  with  the  importance  of  edges  in  medical  images.  If  we  are  to 
include  continuous  images  in  our  space,  then  by  the  Sobolev  Embedding  Theorem 
we  must  use  H®  for  s  >  1.  This  is  also  sufficient  to  guarantee  M®  £  L*®. 

The  requirement  of  this  section,  that  of  a  bounded  restriction  mapping, 
corresponds  to  the  existence  of  a  natural  trace  operator  in  Sobolev  spaces. 

In  the  space  H®(J1)  there  is  a  natural  trace  operator  H®(ft)  H®‘^/^(^).  If 
we  require  7|^  to  be  in  H®(dn),  then  that  too  would  point  to  using  7eH®(li) 
for  s  greater  than  g. 

DIRECT  FORM 

In  this  section  the  Frechet  derivative  of  the  potential  as  a  function  of 
conductivity  is  calculated.  First  the  inverse  operators  are  defined  as 
follows: 

«■*(»)  - 
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is  the  inverse  of  -»  Since  the  Dirichlet  problem 

has  a  unique  solution,  G  is  well  defined  and  the  spaces  have  been  chosen  so 
that  the  inverse  is  bounded.  Similarly 

C":  .  hJcO) 


solve 


V  =  s,  tV„u  =  0,  J 


In  addition  we  will  need  the  mapping 

G®: 

which  solves  the  homogeneous  Neumann  problem 

L/  =  0,  7V  u  =  h,  [  u  =  0 

In  estimating  the  norm  of  which  depends  linearly  on  7,  the  following 
Lemma  will  be  used. 


Proof 


1  ul 


Lemma  1.1  i|L^|  Ig^jj-lll  |7|  |  where  S  is  either  or 

Proof 

liyiliri(«)  =  s“p„6Sch' 

I  L»7Vii.ii-  Lvw.Viil 
'  ^“PyfS  ^  ■■  .V 

=  ®"Pv€S%tV».Vu)/|1w|||,i 

i  IItII  I|u||„i 

We  are  now  ready  to  prove  the  following  theorem: 

Theorem  If  L^u  0  and  l'^+^^(u+^u)  =  0  and  the  Dirichlet  or  Neumann  data 
for  u  and  u  +  ^u  agree  then 

I'^7«  +  yu  =  0(1 1^711^) 

For  the  Neumann  constant  case  we  also  have 
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+  'fV^Sn  =  o(ll^rll^). 

The  forward  napping  defined  by  F(7)  =  n  (which  solves  a  given  Neumann  or 
Dirichlet  problem)  is  C“.  For  the  Dirichlet  case 

r'(T)«T  = 

and  the  Neumann  case 

F'(r)<r  =  -  o5«7V 

Proof 

In  all  cases  we  have 

=  L^tu  +  L^du  +  L^^du  =  0 

First  consider  the  Dirichlet  conditions  constant  case  =  0. 

Applying  G  =  G®  gives 

=  -Clj^u 

Formally 

Sm  =  ^E^(-GL^^)ju 

which  converges  for  ||GL^^||  <  1.  Using  Lemma  1.1 

ll%ll  <  IIGIl-ll^rll 

so  convergence  is  achieved  by  requiring 

ll«7ll 

IM' ' 

(the  norm  used  for  linear  operators  being  the  standard  linear  operator  norm) . 
This  constitutes  a  Taylor  series  for  F,  as  is  linear  in  ^7  the  term  is 
homogeneous  of  order  j  in  Sy.  In  particular  we  have 

r'(i')<T  =  -«%^u 
and 

=  odl^rll^) 

as  claimed.  It  can  readily  be  seen  that  F  is  C®  as  the  higher  derivatives 
can  be  extracted  from  the  Taylor  series. 
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la  tke  ieaaaaB  caasuat.  caae  wt  aill  treat  tJbe  taa  caapoaeats  of 
■•(B)  =  |[^(B)  aefaratdy  ariti^s  T  *  ♦  T^-  'Was  t«a  partial 

4kri«atiaea  caa  be  calcalatd  jepamffly-  la  aMitiaa  to  tbe  epeatiaa 

♦  l/a 

lie  have  its  baaaiary  epaivalcat 


To  calcalate  ve  arrow  =  0. 

The  hoaadary  caoilitiao  aav  refaces  to  if^fa  =  0  aad  the  proof  proceeds  as 

before  with  C  =  <>|[  aad  we  hare 
T 

aad  this 

IF 

la  the  other  haad  IF/ly^  caa  he  calcalated  assaaiag  that  7  =  0  ia  fiP. 
This  leads  to 


L  fa  =  0 

T 

and 


=  -  V»*  -  V.*- 

D 

The  proof  proceeds  in  a  siailar  way  to  the  interior  case.  Applying 

a 

=  C*(-^T^T„n-i7l  Vjjfa) 

rearranging 

The  series 

fm  -  S  (-(^<T,V^" 

converges  for  ||C*«rjV„||  <  1  vhick  is  ensnred  by  Hiyjll  <  l/dlcJlI  IIVjID- 
Thus  we  have  the  derivative 
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aai  tlM!  renlt 


Froa  the  dhaia  rale 

fMh  =  -«^^.  -  cJiT/,.- 

The  faactioa  caa  he  seea  to  be  CT  as  the  partial  derivatives  exist  of  all 
order. 

□ 

la  the  Dirichlet  coastaat  case  F  is  eqaal  to  its  Taylor  expansion  in  a 
nei^boarhood  of  7  aad  is  therefore  aa  aaalytic  napping  (this  was  pointed  out 
Iqr  Calderon  asiag  7=1).  (In  the  Netouan  conditions  constant  case  the  proof 
gives  the  slightly  weaker  result  that  F  is  an  analytic  in  each  coaponent 
separately).  This  indicates  that  the  forward  napping  could  hardly  be  better 
behaved;  however  it  is  the  inverse  of  this  napping  which  is  required  for  EIT 
reconstruction  and  as  we  shall  see  that  is  not  nearly  as  nice. 

nmillETATIM  AS  A  SfUCE 

The  essential  result  of  section  3.4  is  that  to  first  order  we  can  assune 

V.jVSu  =  -V.Sj^n. 

One  interpretation  of  this  is  that  the  perturbed  field  du  is  'caused'  by  a 
distribution  of  current  sources  -V.d7Vu.  Equivalently  one  could  say  that 
adding  a  source  field  s  =  V.dyVu  would  cancel  the  effect  of  changing  the 
conductivity  by  6j. 

An  interesting  case  to  consider  is  to  take  67  =  the  dirac  delta 
distribution  at  a  point  p.  Ve  will  assune  that  u  e  C^(ft)  and  that  p  is  not  a 

critical  point,  Vu(p)  ^  0.  For  sinplicity  we  will  take  7=1.  The  source 
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ten  is  now  -  which  is  a  dipole  with  dipole  nonent  | Vo (p)| oriented  in 

the  direction  of  the  current  vector  at  p.  The  perturbation  du  to  first  order 
is  then  the  electric  field  associated  with  the  dipolee.  The  function 
is  the  point  response  of  the  systea  (up  to  first  order) ,  in  optics  this  would 
be  called  the  point  spread  function.  In  contrast  to  ideal  optical  systems 
the  response  is  position  dependent  falling  off  dramatically  as  p  gets  further 
from  the  boundary.  The  field  from  a  dipole  is  asymptotically  (2/t  cos  ^)/r  . 
Here  n  is  the  dipole  moment,  r  the  distance  from  the  dipole  and  6  the  angle 
relative  to  the  dipole  orientation.  In  this  case  /i  =  |Vu|  which  is  at  best 
constant  and  typically  decreases  away  from  Xl.  Hence  we  find 
1  lH^/^(dQ  ■  )  ''**®*^®  P  =  dist(p,dft).  This  means  that  the  ability 

of  an  EIT  system  to  detect  a  small  object  of  high  conductivity  contrast  will 
fall  off  at  best  proportionately  to  the  inverse  square  of  distance  from  the 
boundary. 

INTEGRAL  FORM 

It  is  useful  to  reformulate  the  linearised  problem  in  an  integral  form.  In 
Chapter  4  the  finite  element  method  will  be  used  to  represent  the  electric 
potential  and  to  solve  the  forward  problem.  In  the  finite  element  method 
differential  equations  are  formulated  as  variational  problems,  this  is 
equivalent  to  the  weak  form  of  the  differential  equation.  Since  this  is 
essentially  an  integral  formulation  it  will  be  advantageous  to  express  the 
linearised  conductivity- to- voltage  mappin  as  an  integral  operator. 

1 

First  notice  that  for  any  U  with  L^U  =  0  and  any  V  (U  and  V  in  H  (fi)  say) 

V.(U7VV)  =  7VU,VV 

and  so  using  the  divergence  theorem 
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Froa  above 


f  V7V„U  =  [  rW.VU. 

JaQ  “  Jfl 


=  o(jl^rlr) 


where  L^u  =  (u+tfu)  =  0  and  =  (r+^7)Vjjn,  Choose  any  v  with 


b^w  =  0  then 
7 


V.w(']'VSn+SjVu)  =  7Vw.V5u+67Vw.Vn+o(|  |^7|  I  ). 
Applying  the  divergence  theoren 

J^w(7Vn^u+^7VnU)  =  J  7?w.V^u+^cVw.Vu+o(l  |^7|  |^) 


Since  L^w  =  0  we  know  that 


I  ^U7V„w  =  [  7V5uVw 
JdQ  ®  h 

also,  from  the  boundary  conditions  +  SyVnV^  =  o(l  l^7l  I  )>  result  is 

J^^U7Vj^w  =  -J  ^7Vw.Vu  +  o(|  |^7|  |^). 

Let  us  now  consider  the  implications  of  this  formula  for  the 
reconstruction  process.  One  has  some  initial  estimate  of  the  conductivity  7 
and  wishes  to  correct  this  using  the  best  linear  approximation.  Some  known 
current  patterns  are  applied  to  the  surface  of  the  body  Xl.  Measurements 
of  voltage  u.  are  made  between  various  electrodes.  Since  measurement  is  an 
averaging  process  over  the  electrode  we  will  assume  that  the  measurements  are 
of  the  form 

where  Wj^  is  characteristic  of  the  geometry  and  electrical  properties  of  the 
electrode  pair  k.  We  have  an  a  priori  model  of  the  body  with  conductivity  7 
which  we  compare  with  the  real  body  which  has  conductivity  7  +  ^7.  The 
discrepancy  between  the  two  is  measured  by  the  data 

CTi_k  =  f 
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We  thea  solve  the  Neunann  problea  L^w  =  0  subject  to  Vj^Wj^  ”  ^ 

linear  approxiaation  to  the  conductivity  error  S'!  we  solve  the  system  of 
linear  equations 

*™i,k  = 

In  this  formulation  the  Neumann  conditions  were  kept  constant.  This  is 
the  most  useful  formulation  for  impedance  measurement  for  both  theoretic  and 
practical  reasons.  For  completeness  we  must  compare  this  to  the  problem 
investigated  by  Calderon  in  which  the  Dirichlet  conditions  were  fixed  and  a 
difference  of  Neumann  conditions  (that  is  boundary  current  densities) 
measured.  In  this  case  the  boundary  conditions  are  j  =  -^V^^u  and 

assuming  now  that  =  0.  This  leads  to  the  result 

[  -  f  ^tVw.Vu  +  0(11^711^) 

which  does  not  have  the  curious  minus  sign 

Not  a  radon  transform 

It  has  been  assumed  by  a  number  of  authors  (Schaumberg  [11] ,  Barber  and  Brown 
[12],  Vogelius  and  Santoza  [13]  that  the  mapping  F:^7  ^  ^u|^  can  be 
approximated  by  a  Generalised  Radon  Transform  (GRT) ,  in  the  sense  of  Quinto 
[14].  Schaumbert  assumed  that  the  fibres  were  the  current  streamlines 
whereas  Barber  and  Brown  and  Vogelius  and  Santoza  assumed  them  to  be  the 
equipotential  lines  (in  the  case  ficR  ).  Such  GRTs  are  linear  maps  so  it  is 
reasonable  to  compare  them  with  the  Frechet  derivative  F(7)  which  is  in  a 
specific  sense  the  best  linear  approximation. 

Consider  as  before  the  unit  disc.  Take  J  =  cos  7  =  1.  The  solution 
of  the  Neumann  problem  for  Laplace's  equation  is  clearly  n{r,9)  =  r  cos  ^  or 
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=  ( 

'  it 


equivalent  u(x,y)  =  x.  If  the  perturbation  ^7  is  circularly  symmetric  then 

R 

cos  0  will  still  be  an  eigen  funtion  of  1  +  6j.  For  simplicity  take 

lxl>/) 

|xl<p  • 

then  ^  “  A^cos  6  and  Sn  =  5j'(0)cos  9  +  o(|5/?p).  Thus  ^u  has 

support  on  3ft  -  {T/2.3r/2}.  If  F'(l)  were  a  GRT  then  the  support  of  3u  would 
be  contained  in  the  'shadow'  of  the  set  S  =  {x:|x|</?},  that  is  all  the  3  6  3ft 
such  that  <5^  n  S  O  where  <5^  is  the  fibre  through  9.  It  is  clear  that  this 
is  not  the  case  in  this  example  either  for  the  current  streamlines 
=  {(x,y):y=cos  3}  nor  for  equipotential  lines  ^  =  {(x,y):x=cos  3}. 


Reconstruction  algorithms 

Some  of  the  reconstruction  algorithms  suggested  by  this  work  have  been 
implemented  and  numerical  results  are  reported  in  Breckon  and  Pidcock  {15] . 

We  will  only  give  brief  details  here.  The  general  procedure  is  as  follows. 
Make  an  initial  guess  7^^^  to  7 
Repeat 

Choose  current  patterns  > ^2^”^ ’ "  ’ ’ ^m^"^ ’  possible  using  an 
optimal  procedure  and  make  measurements  3V^j. 

Find  a  least  squares  solution  to 
7Vwj^.Vu.  for  all  k,i 

Update  =  7(**)  +  6j 

Until  3V.j^  <  G 

The  linearised  inverse  problem  is  extremely  ill  posed  as  can  be  seen  from 
the  singular  values  given  in  [16] .  Consequently  regularisation  must  be  used 
in  solving  the  linear  system.  If  only  one  iteration  were  used,  the  small 
number  of  useful  data  values  given  by  the  optimal  methods  (i.e.  only  those 
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which  can  distinguish  7^“^  from  7)  would  be  a  severe  problem.  However  in 
both  of  the  adaptive  methods  discussed  new  currents  are  found  at  each 
iteration  to  distinguish  best  between  what  is  really  there  and  the  latest 
guess.  The  map  will  be  completely  different  at  each  iteration 

giving  rise  to  completely  different  current  patterns  and  voltage  data.  While 
the  data  set  measured  at  each  iteration  is  typically  smaller  than  that  used  by 
non- adaptive  methods,  the  information  content  of  the  data  is  maximised  at  each 
step. 
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S.R.  CLOUDE 

Polarisation  in  electromagnetic  inverse 
problems 

ABSTRACT 

With  recent  advances  in  measurement  technology,  full  vector  information  is  now 
available  from  electromagnetic  scattering  experiments.  This  paper  presents  a 
review  of  techniques  developed  for  the  inclusion  of  this  vector  information  in 
inverse  scattering  theory. 

The  magnetic  field  integral  equation  is  used  to  develop  three  inverse 
scattering  models  and  to  demonstrate  the  importance  of  full  polarisation 
information  for  3D  reconstruction.  We  also  consider  two  other  inversion 
techniques,  vector  diffraction  tomography  and  inverse  boundary  conditions,  and 
discuss  the  polarisation  aspects  of  each. 

1.  INTRODUCTION 

Very  often  the  vector  nature  of  electromagnetic  waves  is  ignored  in  direct 
scattering  problems.  Solutions  are  then  formed  for  a  complex  scalar 
wavefield  satisfying  the  scalar  Helmholtz  equation.  This  approximation  is 
justified  only  when  2- dimensional  problems  are  considered  (in  which  case  we 
can  treat  TE  and  TM  waves  separately)  or  when  paraxial  solutions  are  adequate 
and  crosspolarisation  is  not  of  interest.  A  further  advantage  of  adopting  a 
scalar  approach  is  the  unity  it  provides  with  acoustic  scattering  theory  (see 
Jones  1986). 

For  3-D  inverse  problems  however,  this  scalar  approximation  is  neither 
justified  nor  desirable.  Measurement  techniques  have  been  developed  (mainly 
as  a  result  of  interest  in  frequency  re-use  in  communications)  to  accurately 


measure  vector  field  quantities  (see  the  review  by  Cox  1981).  Here,  we 
address  the  problem  of  how  such  information  can  be  used  to  help  in  the 
inversion  process.  This  provides  a  great  challenge  to  inverse  scattering 
theory;  a  complete  inversion  of  the  vector  problem  still  evades  us  (see 
Langenberg  1989)  but  many  of  the  concepts  required  have  been  identified. 

They  involve  advanced  ideas  from  vector  scattering  theory,  as  well  as  such 
unlikely  elements  as  Lie  algebra  and  group  theory.  While  the  subject  is 
still  in  the  early  stages  of  development,  we  present  here  a  selective  review 
of  techniques  of  use  in  the  vector  inverse  problem.  For  the  sake  of  brevity, 
we  concentrate  only  on  those  features  which  relate  directly  to  polarisation  of 
the  scattered  field,  leaving  other  details  to  the  references.  This  may  be 
considered  an  update  of  a  similar  review  published  in  1981  by  Boerner. 

With  this  in  mind  we  develop  the  paper  in  three  main  stages:  first  we 
review  the  matrix  algebra  used  to  describe  polarisation  effects  and  consider 
general  symmetry  properties  of  the  scattered  field  which  follow  from  generic 
properties  of  the  scatterer.  These  are  important  because  they  may  be  used  to 
impose  global  or  local  symmetry  on  the  reconstruction.  We  then  consider  a 
set  of  inversion  techniques  based  on  the  magnetic  field  integral  equation 
(MFIE),  with  the  high  frequency  physical  optics  and  extended  physical  optics 
theories  as  special  limiting  cases.  In  particular,  we  review  the  important 
Kennaugh- Cosgriff  inversion  formula  and  show  how  we  may  add  a  polarisation 
correction  term  (first  derived  by  Bennet)  to  yield  information  on  specular 
point  curvature.  We  then  describe  an  exact  inversion  method  based  on  the 
MFIE  and  discuss  polarisation  aspects  of  two  other  inversion  schemes;  vector 
diffraction  tomography  using  the  dyadic  Greens  function  (Langenberg  1989)  and 
the  inverse  boundary  condition  method  as  developed  by  Imbriale  and  Mittra 
(1970) . 
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2  MAXm  FOUIAIISM 

Iftider  a  far  field  assumption,  we  may  write  the  vector  scattered  field  as  a 
multiple  convolution  between  the  incident  field  and  object  impulse  response  as 

®M  =  %  *  ®N  0  ^  <  1  (1) 

By  a  Fourier  transform  we  obtain  a  2x2  complex  scattering  matrix  [S]  at 
angular  frequency  u  as 

Ss=ra.Ei  (2) 

where  and  Ej  are  spinor  quantities  representing  the  scattered  and  incident 
electric  fields  respectively.  The  spinor  nature  of  the  electric  field 
follows  from  its  transformation  properties  under  a  change  of  (complex) 
orthogonal  base  states  (see  Cloude  1986).  We  may  write  such  a  transformation 
as 

S'  =  [UJ-S  (3) 

where  [U2]  is  a  2x2  unitary  matrix  with  unit  determinant.  We  may  interpret 
this  geometrically  as  a  rotation  in  a  real  three-dimensional  space  by  writing 
[U2]  as 

[U2]  =  cos  ^  (Tq  -  i  sin  $  {<rj^  cos  a  +  try  cos  +  o-gcos  7) 

=  exp(-i<?[£.n])  (4) 

representing  a  rotation  of  26  about  an  axis  specified  by  angles  0,  /?,  7.  The 
set  £  =  <^re  the  Pauli  matrices.  From  spinor  algebra  we  know 

that  we  may  associate  with  (3)  a  spin  matrix  or  quaternion 

[X]  =  +  yiTy  +  (5) 

where  the  real  3  vector  r  =  (x,y,z)  maintains  unit  modulus  under  a 
transformation 

[X]  =  [U2][X][U2]"  (6) 

or 

r'  =  [OgJ.r  (7) 
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where  [Og]  is  a  3*3  real  orthogonal  matrix.  The  locus  formed  under  unit 
modulus  is  termed  the  Poincare  sphere  (see  Born  and  Volf  1965),  which  sits  not 
in  physical  space,  but  in  an  abstract  polarisation  space  of  the  spinor  E. 

Many  measurement  systems  operate  by  measuring  the  vector  r  rather  than  the 
complex  spinor  E  (Cloude  1989).  The  two  are  related  by  the  spin  matrix  X  as 

X  =  EjEj*  -  EyEy* 
y  =  Re(EjEy*) 
z  =  Im(EyEy*) 

If  we  further  define  m  =  Ej^E^*  +  EyEy*  then  g  =  (m,r)  is  called  the  Stokes 
vector  of  the  wave  spinor  and  [X]  =  m<rQ  +  v.e  is  called  the  wave  coherency 
matrix. 

In  scattering  theory,  we  generalise  the  above  by  considering  (2)  as  a 
spinor  transformation 

[S]  =det([S])[m[H] 

=  d.exp[(-ife+Aii).£]  (8) 

where  [E]  is  2x2  unitary  and  [H] ,  2x2  hermitian  (we  shall  assume  for  the 
moment  that  det([S])  is  nonzero,  it  becomes  zero  only  for  degenerate 
scattering  systems  such  as  a  linear  dipole  or  helix).  This  transformation 
corresponds  to  a  combination  of  boost  and  rotation  of  the  Stokes  vector  and 
has  8  degrees  of  freedom;  2  direction  vectors  n  and  p  lying  in  the  space  of 
the  vector  r^  and  2  angles  ^  and  A.  The  direction  vectors  n  and  p  figure 
prominently  in  the  theory  of  null  polarisations  as  developed  by  Kennaugh  (see 
Kennaugh  1952  and  Boerner  ''81).  He  showed  that  for  backscatter,  [S]  has  two 
orthogonal  eigenstates  and  two  copolar  null  states  (where  the  incident  wave  is 
scattered  into  an  orthogonal  state).  Since  these  states  display  special 
symmetries  of  the  scattering  volume  they  have  been  proposed  as  important 
features  in  the  inverse  problem  (Boerner  1981),  but  to  date  no  successful  use 


has  been  made  of  these  null  states,  mainly  because  of  the  complexity  of  their 
dynamics. 

By  using  the  Stokes  vector  g  instead  of  the  spinor  E  we  may  express  (2) 
and  (8)  in  the  form 

=  M  El  (9) 

where 

[K]  =  I  Tr(£  S  £  S"^)  =  ®  S'^A  (10) 

[M]  is  called  the  Mueller  or  Stokes  reflection  matrix,  Tr  stands  for  the  trace 

of  the  matrix  and  A  is  a  4x4  unitary  matrix  of  the  form 

1  1  0  O' 

0  0  1  -i 

0  0  1  i 

1  -1  0  0. 

The  matrix  [M]  is  very  important  in  vector  scattering  theory;  the  key  problem 
in  inverse  scattering  is  to  relate  the  elements  of  [M]  to  shape  and  material 
properties  of  the  scattering  volume.  Ve  also  note  that  the  Stokes  vector  g 
may  also  account  for  partially  polarised  waves  where  m  >  r.r.  Such  states 
lie  inside  the  Poincare  sphere  (with  random  polarisation  at  the  origin). 

This  extra  degree  of  freedom  in  g  means  that  [M]  is  more  general  than  [S] ,  it 
contains  all  the  relative  phase  information  of  [S] ,  together  with  information 
on  correlations  between  various  elementary  scattering  mechanisms.  To  see 
this  we  outline  an  alternative  matrix  formalism  based  on  a  4x4  hermitian 
covariance  or  coherency  matrix  (see  Gloude  1986,  1989).  The  covariance 
matrix  [S]  is  related  to  m--,  the  elements  of  [M] ,  by  16  Dirac  matrices, 
formed  from  direct  products  of  the  Pauli  matrices  as 

[E]  =  mj.  (T.  ©(Tj  (11) 

The  scattering  coherency  matrix  [T]  is  based  on  complexification  of  the  spin 
matrix  X  and  is  related  to  [M]  by  a  set  of  16  traceless  hermitian  matrices 
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as 

m  =  •ij  «,i.j  («) 

viiere 

^4i*j  = 

The  Matrices  [£]  aad  [1]  traasfoni  mder  a  SMilariti’  traasfonatiM  iavolvii' 
a  4x4  naitary  Matrix  [lE^] ,  represeati^  a  ckai^e  i>  ceaplex  Matrix  basis  for 
tbe  expaasiM  of  [S],  i.e.  if  [S]  =  k.r  tiuni 

r  =  i%n  (14) 

and 

[T-]  =  [tjm  (>!,]*  0-5) 

where  [T]  =  k-k^-  For  exaqile.  we  aay  deteniiae  the  fora  of  [S]  nadcr 

baclscatccr  when  the  object  has  X  fold  spaetry  in  a  plane  perpend icalar  to 

the  line  of  sight.  FroM  reciprocity  it  follows  that  [S]  is  synoctric.  i.e. 

Sjiv  =  Sy||,  [N]  is  SYMnctric  and  T  and  £  are  3x3  heraitian.  Ihider  a  plane 

rotation  t.  k  becoMCS  k'  =  [11^]  .k  where  [0^]  is  given  by 

i  0  0  0 

0  cos2r  -sin2r  0 

0  sin2r  cos2r  0 

pool 

Since  =  0  for  backscatter  it  follows  that  for  objects  with  N  fold  syaactry 
(and  N  >3)  [S]  Must  have  the  fora  of  a  coaplcx  scalar,  i.e. 

P'l  =  kQ  iTq  (16) 

This  is  a  sinple  exanple  of  what  can  be  done  by  considering  synactry 
properties  of  Matrix  descriptors  (see  Van  de  Hulst  1981  for  further  discussion 
of  symmetry  constraints  on  [M]).  We  can  generalise  such  argunents  to 
arbitrary  unitary  transformations  of  k,  when  we  need  to  involve  Lie  algebra  to 
ascertain  invariant  features  under  elementary  rotations.  Further  details  may 
be  found  in  the  references  (Cloude  1986,  1989).  We  now  turn  to  scattering 
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theory  ia  sb  attempt  to  relate  tlicse  natrix  observables  to  shape  and  naterial 
properties  of  the  scatterin'  volnae. 


3.  ELECnnUCKEriC  IRVEISE  SCATTaiNC 

la  this  section  ve  consider  the  problca  of  reconstruction  of  a 
three- dimensional  perfectly  conducting  body  nsing  the  magnetic  field  integral 
equation,  the  vave  equation  vith  inverse  boundary  conditions  and  vector 
diffraction  tomography.  In  particular,  vc  show  how  full  polarisation 
infomation  is  needed  for  accurate  reconstruction  when  the  object  has  unequal 
principal  cunatures,  and  show  how  such  curvature  inforoation  is  contained 
within  the  elements  of  the  Mueller  matrix  [N] . 

Vc  begin  with  Maxwells  equations  in  the  time  donain, 

V  X  E(r,t)  =  -p  S}l(r,t)/St 

V  X  II(r,t)  =  €  6E(r,t)/a  +  J(r,t) 

V.H(r,t)  =  0 

V.E(r,t)  =  p(r,t)/e 

where 

II(r,t)  =  If II  V  X  A(r,t) 
and 

A(r,t)  =  /{/4x  f  J(r',t-T)/R  dV 
Jy 

Tie 

By  defining  the  total  magnetic  field  as  11  =  H  +  H  and  imposing  the 

T 

boundary  condition  n  x  II  =  Js,  we  obtain  an  integral  equation  (the  MFIE)  for 
surface  current  Js  as  (see  Poggio  and  Miller  1973) 

Js(r,t)  =  2n  X  +  l/2x  n  x  [  L(Js)  x  R  ds  (17) 

where  L  is  a  differential  operator 


1  6  1 


Once  the  current  is  Known,  the  scattered  far  field  may  be  determined  from 

63s 

X  r  ds  (18) 


IF 


Ve  may  write  the  MFIE  (17)  as  the  sum  of  three  contributions: 

J(r,t)  =  2n  X  +  pJ(r,t)  +  S  (19) 

where  the  first  term  represents  the  direct  influence  of  the  incident  field, 
the  second  represents  a  self  patch  contribution  to  the  current  and  the  third, 
S,  is  an  integral  over  other  patches  with  currents  evaluated  at  earlier  times. 
The  factor  p  has  been  shown  by  Mieras  and  Bennet  (1982)  to  be  (Marx  (1985) 
derived  a  more  general  expression  for  the  self  patch  term  which  reduces  to  the 
Bennet/Mieras  result  for  constant  current  across  the  patch) 

1  Im. 

o  =  Ky)  (20) 

where  Ky  and  Ky  are  principal  curvatures  and  AA  is  the  area  of  the  patch. 

Ve  may  now  derive  inverse  scattering  identities  by  considering  various 
approximations  to  the  surface  current.  In  the  simplest  case  we  ignore  the 
self  patch  term  and  the  integral  S  and  consider  only  the  forcing  term.  If  we 
further  assume  the  current  to  be  zero  on  the  shadow  side  of  the  object  we 
obtain  the  well-known  physical  optics  assumption.  The  expression  for  the 
scattered  far  field  (18)  then  becomes 

ng(r,t)  =  A  6/6t  (n  x  Ijp  x  r  ds  (21) 

If  we  assume  an  incident  impulse  plane  wave  with  E  polarised  in  the  x 
direction  and  travelling  in  the  positive  z  direction  then  the  vector  triple 
product  reduces  to  a  simple  area  projection  funtion  which,  when  Fourier 
transformed  and  integrated  by  parts,  yields  the  following  expression  for  the 
impulse  response  (see  Kennaugh  and  Moffat  1965) 
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where  A(t  )  is  the  projected  area  funtion  of  the  target,  i.e.  the  cross 
sectional  area  of  the  target  in  a  plane  perpendicular  to  the  z  axis, 
delineated  by  a  plane  moving  with  a  speed  of  one  half  light  velocity.  For 
example,  for  axial  incidence  on  a  rotationally  symmetric  target  k  =  rp  ,  where 
p  is  the  target  contour  function.  This  important  result  was  first  derived  by 
Kennaugh  and  Moffat  in  1965  and  has  been  used  by  several  authors  as  the  basis 
for  electromagnetic  imaging  (see  Young  1976,  Shubert  1977).  Note  that  while 
the  early  time  impulse  response  is  well  approximated  by  this  technique,  late 
time  ringing  due  to  damped  current  flow  on  the  object  is  ignored  (it  is 
contained  in  the  integral  S  which  we  chose  to  ignore).  This  technique  has 
been  extended  to  account  for  late  time  damping  using  GTD,  moments  of  the 
impulse  response  and  phenomenological  damping  terms  (Kennaugh  and  Moffat 
1965) . 

The  technique  known  as  ramp  response  imaging  (Young  1976)  is  based  on 
integrating  the  impulse  response  twice  to  yield  a  direct  relationship  between 
a  measured  waveform  and  the  area  function  A(t„).  This  has  been  used  to 
derive  information  on  object  cross-section  as  a  function  of  distance  along  the 
line  of  sight  as  well  as  object  length  and  total  volume  (from  the  integral  of 
the  ramp  response).  By  using  multiple  looks,  limiting  contours  have  been 
found  for  object  shapes  such  as  cubes  and  concentric  cylinders.  Note  that 
this  time  domain  result  is  related  to  the  well-known  Bojarski  identity  by  a 
Fourier- Radon  transform  (as  shown  by  Boerner  1978). 

For  our  puposes  we  note  that  under  physical  optics  there  is  no 
depolarisation  of  the  backscattered  signal.  However,  this  does  not  mean  that 


there  is  zero  depolarisation  in  the  high  frequency  limit  of  the  MFIE,  as  we 
not  only  ignored  the  integral  S  but  also  the  self  patch  term  in  deriving  this 
result.  Since  the  latter  is  proportional  to  the  radius  of  the  self  patch,  we 
might  expect  it  to  be  of  only  second  order  importance  (very  often  the  self 
patch  term  is  ignored  in  numerical  calculations  using  the  MFIE).  However, 
this  is  not  the  case  for  general  scattering  bodies.  When  determining  the 
scattering  matrix  for  the  object,  the  self  patch  term  makes  an  important 
contribution  to  the  phase  difference  between  copolar  terms,  as  well  as 
generating  finite  cross  polar  terms.  In  fact,  the  self  patch  term  is  only 
zero  for  surfaces  which  have  equal  radii  of  curvature,  such  as  a  flat  plate  or 
sphere. 

If  we  assume  the  current  over  the  self  patch  to  be  of  constant  magnitude 
and  equal  to  the  physical  optics  current  then  we  may  write  a  better 
approximation  to  the  surface  current  as  (Mieras  and  Bennet  1982,  Bennet  1978) 

J  =  2nxf  *€  (K„-Kv)(Jpo^u-  JpoV»/4  (23) 

where  we  have  assumed  a  circular  patch  of  radius  6  =  ct-z  and  the  unit  vectors 
u  and  V  are  aligned  with  the  principal  coordinate  axes  of  the  surface  at  the 
point  of  interest.  When  we  integrate  this  expression  to  obtain  the  scattered 
field,  the  e  factor  yields  a  term  proportional  to  the  first  instead  of  second 
derivative  of  the  area  function.  The  final  result  for  the  impulse  response 
is  of  the  form  (see  Foo  1984,  Chaudhuri  1986,  1977) 

1  A  .  ^  %%  M 

'■  =  2?  ^  TTJt  ((%■«)"■  (Sii-xW  (2-*) 

where  we  have  assumed  a  horizontally  polarised  incident  magnetic  field  with 
unit  vector  ajj.  We  can  now  determine  the  elements  of  the  scattering  matrix 
[S]  as 
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1  S^k  6k 


(25) 


=  2t  ^^2  4t  H 

hh  Sk 
\\  ^  “IT  'Ei 

1  6k 

^vv "  ^  ^  “4r  ^ 

where  a  is  the  angle  between  the  wave  coordinates  and  principal  coordinates 

for  the  specular  point.  If  we  expand  this  matrix  in  terms  of  the  Pauli 

matrices  and  then  perform  a  Fourier  transform  with  respect  to  k^,  the  free 
space  wavenumber,  we  obtain  a  complex  vector  k  of  the  form 

k  =  (kQ^A(kQ)/2x,  A  ikQA(kQ)cos  2a,  A  ikQA(kQ)sin  2a,  0)  (28) 

where  A  is  the  difference  in  curvature  at  the  specular  point  and  A(kQ)  is  the 

value  of  the  Fourier  transform  of  the  area  function  at  spatial  wavenumber  kQ. 
Note  that  A  is  real,  so  the  second  two  elements  are  in  phase  quadrature  with 
the  first.  This  implies  that  Sjjjj  and  Syy  have  equal  magnitude.  Their  phase 
difference,  however,  is  directly  related  to  the  factor  A.  This  means  that  we 
can  estimate  curvature  difference  at  the  specular  point  by  measurement  of  the 
phase  difference  between  copolar  terms  (measured  at  a  wavelength  much  shorter 
than  body  dimensions,  see  Chaudhuri  1986,  Boerner  1981).  Note  we  also  have 
information  on  the  size  of  the  object  through  the  factor  A(kQ)  and  on  local 
orientation  through  the  angle  a. 

Ve  can  avoid  the  need  for  coherent  measurements  (accurate  phase 
measurements  are  very  difficult  to  achieve  in  practice)  by  calculating  the 
coherency  and  Mueller  matrices.  The  three  pieces  of  information  A(kQ),  A  and 
a  are  then  encoded  in  the  diagonal  elements  of  [M]  as 

=  4Ar([M])  =  4r^(mQQ+mjj+m22+iii33)  (29) 
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(30) 


=  2k/(«(,„  -  n,33)/Tr([M]) 


tan  2a  = 


W"ll'"'22""'33 


The  12  off-diagonal  terms  of  the  Mueller  matrix  relate  to  correlations  between 
curvature  difference  A,  amplitude  weighting  A(kQ)  and  orientation  a  and  hence 
for  single  point  scattering  provide  no  extra  target  information.  For  complex 
objects  however  (e.g.  rough  surfaces)  these  additional  terms  are  useful  for 
assessing  variation  of  surface  topology. 

The  next  logical  stage  is  to  include  the  whole  MFIE  in  the  determination 
of  surface  current.  Unfortunately  we  cannot  then  obtain  analytical  solutions 
but  must  resort  instead  to  numerical  techniques.  Nonetheless,  Bennet  (1981) 
has  used  the  full  time  domain  MFIE  for  object  reconstruction  using  two 
numerical  techniques,  the  first  based  on  iteration  using  successive  estimates 
of  the  surface  current  and  measured  ramp  response  and  a  second  based  on  a 
straightforward  extension  of  the  time  stepping  procedure  used  in  direct 
scattering  implementations  of  the  MFIE.  However,  examples  are  only  worked 
for  axial  incidence  on  bodies  of  revolution  where  the  area  function  may  be 
parameterised  in  terms  of  a  contour  function  p  and  depolarisation  does  not 
occur.  Although  extension  to  non- symmetrical  bodies  was  attempted,  this  area 
still  remains  to  be  fully  investigated. 

In  summary,  we  have  seen  that  we  can  obtain  three  inverse  schemes  based  on 
the  MFIE,  the  first  yielding  the  physical  optics  Kennaugh- Cosgrif f  formula 
which  is  a  time  domain  version  of  the  Bojarski  identity.  Ve  then  saw  that  by 
including  the  self  patch  term  we  obtained  a  polarisation  dependent  correction 
to  the  inversion  formula  whereby,  in  the  high  frequency  limit,  the  phase 
difference  between  copolar  terms  is  simply  related  to  curvature  difference  and 
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depolarisation  is  caused  by  misalignment  of  wave  and  surface  coordinates. 
Finally  we  saw  how  the  full  MFIE  may  be  used  for  object  reconstruction  by 
using  numerical  techniques  akin  to  the  time  stepping  procedures  well  known 
from  direct  scattering  applications  of  the  MFIE. 

3.1  Inverse  boundary  conditions  and  vector  diffraction  tomography 
As  an  alternative  to  the  MFIE  formulation  of  electromagnetic  inverse 
scattering,  Imbriale  and  Mittra  (1970)  devised  a  technique  based  on  an  inverse 
boundary  condition  for  the  wave  equation.  Weston  and  Boerner  (1969)  showed 
that  this  condition  (namely  that  the  total  tangential  electric  field  must  be 
zero  at  a  surface  of  a  perfect  conductor)  is  sufficient  to  reconstruct  the 
object. 

In  this  technique,  complex  measured  data  is  required  over  an  enclosing 
sphere  at  one  frequency.  This  data  is  then  Fourier  transformed  to  obtain 
weighting  coefficients  in  a  Fourier  expansion  in  terms  of  the  angular  variable 
0.  These  coefficients  are  then  used  together  with  the  appropriate  Greens 
function  (a  Hankel  function  for  the  2-D  problemr  considered)  to  seai'ch  for  a 
point  at  which  the  total  field  (i.e.  incident  plus  scattered)  is  zero  (or  some 
minimum  if  limited  aspect  data  are  obtained).  The  technique  effectively 
involves  analytic  continuation  of  the  field  up  to  the  surface  of  the  scatterer 
and  can  be  modified  for  concave  bodies  and  multiple  scattering  (see  Mittra 
1973) . 

Unfortunately,  although  the  technique  is  perfectly  general  and  can  in 
principal  be  used  for  the  three-dimensional  case,  only  two-dimensional 
examples  have  been  published.  Since  2-D  problems  can  always  be  decoupled 
into  TE  and  TM  waves,  the  significance  of  polarisation  for  the  3-D  case  has 
not  been  clearly  developed.  Note  that  this  technique  is  a  point  by  point 


reconstruction  and  such  tends  to  be  computationally  intensive.  However,  for 
many  applications  only  a  limited  number  of  surface  points  may  be  required  and, 
with  the  widespread  availability  of  super  computer  power,  this  technique  may 
become  feasible  for  low  to  intermediate  frequency  applications  (the 
Kennaugh- Gosgrif f  formula  being  more  efficient  for  high  frequency  problems). 
Ahluwalia  and  Boerner  (1973)  considered  a  generalisation  of  this  technique  to 
lossy  dielectric  bodies  and  showed  that,  while  uniqueness  is  no  longer 
assured,  useful  reconstructions  can  still  be  made,  given  limited  a  priori 
target  information. 

The  third  technique  of  current  interest  for  3-D  inverse  scattering  is  that 

based  on  a  vector  extension  of  linearised  diffraction  tomography  (see 

Langenberg  1989).  In  this  technique,  the  vector  Greens  theorem  is  used  to 

obtain  an  electromagnetic  version  of  the  Porter- Bor jarski  integral  equation 

for  a  vector  holographic  field  6  as 

00 

9  =  -2up  jjj  Jc(R,w)Gj(R-R'  ,w)d^R  (32) 

-00 

Vhere  G^  is  the  imaginary  part  of  the  vector  Greens  Function  given  by 

G  =  (I+W)g  (33) 

where  I  is  the  identity  dyadic  and  g  the  scalar  free  space  Greens  function. 

In  order  to  invert  the  Porter- Bojarski  formula  we  must  perform  two  operations; 
an  inversion  of  the  convolution  (as  in  the  scalar  case)  plus  the  added  problem 
of  inversion  of  the  dyadic  G^.  The  former  can  be  handled  by  linearising, 
assuming  physical  optics,  and  integrating  with  respect  to  the  frequency 
variable  w.  The  latter,  however,  is  difficult  because  the  imaginary  part  of 
the  Greens  dyadic  has  zero  determinant  (in  far  field  scattering  it  effectively 
projects  3-D  fields  onto  a  transverse  plane).  The  result  is  that  only 
projections  of  the  solution  can  be  calculated,  leading  to  two  simultaneous 
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scalar  inversion  integrals  for  a  dyadic  V7,  where  7(R)  is  the  object  singular 
function  which  is  related  to  the  characteristic  function  considered  in  the 
traditional  scalar  Bojarski  theory.  The  net  result  is  that  differential 
geometric  manipulations  are  needed  to  obtain  reconstructions  using  this  vector 
theory  (see  Langenberg  1989) .  The  relationship  between  diffraction 
tomography  and  the  extended  physical  optics  theory  of  Bennet  has  yet  to  be 
evaluated. 

A  SUMMARY  AND  CONCLUSIONS 

In  this  paper  we  have  reviewed  techniques  of  use  in  the- formulation  of  vector 
inverse  scattering  problems.  We  have  seen  that  matrix  algebra  may  be  used  to 

provide  a  formalism  suitable  for  measurement  and  analysis  of  polarisation 
problems  and,  in  particular,  that  transformation  properties  of  the  matrix 
descriptors  may  be  used  to  establish  symmetries  of  the  scattering  volume.  We 
have  also  seen  a  clear  relationship  between  the  coherent  scattering  matrix  and 
real  Mueller  matrix,  with  the  latter  providing  information  of  correlation 
properties  in  complicated  scattering  scenarios. 

The  magnetic  field  integral  equation  may  be  used  to  establish  three 
important  inversion  techniques  by  using  three  approximations  for  the  surface 
current.  The  importance  of  polarisation  is  clearly  demonstrated  by 
considering  the  effect  of  the  self  patch  term  in  the  MFIE.  It  leads  to  a 
simple  but  important  relationship  between  the  copolar  phase  difference  and 
surface  curvature. 

Other  electromagnetic  inverse  techniques  which  show  promise  for  3D 
inversion  are  those  based  on  inverse  boundary  conditions  and  the  point  by 
point  reconstruction  technique  of  Imbriale  and  Mittra  and  those  using 
manipulations  of  dyadic  Greens  functions  in  a  generalisation  of  the  Porter 
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Bojarski  integral  equation. 
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G.  DASSIOS  AND  K.  KIRIAKI 

Size  orientation  and  thickness  identification 
of  an  ellipsoidal  shell 

ABSTRACT 

A  confocal  ellipsoidal  shell  of  unknown  seiniaxes,  orientation  and  thickness  is 
excited  into  secondary  radiation  by  a  plane  harmonic  wave  at  low-frequency. 

The  outer  ellipsoid  forms  a  penetrable  surface  while  the  inner  one  is  a  soft 
ellipsoid.  We  show  that  one  measurement  of  the  leading  low-frequency 
approximation  and  six  specific  measurements  of  the  next  approximation  for  the 
scattering  cross-section  are  enough  to  determine  the  size,  the  orientation  and 
the  shell  thickness  of  the  target. 

1.  INTRODUCTION 

In  a  recent  paper  [4] ,  the  first  author  has  developed  a  simple  algorithm  that 
provides  the  size  and  the  orientation  of  an  unknown  soft  ellipsoid  out  of 
seven  measurements  of  the  forward  scattering  amplitude.  Actually,  only  six 
measurements  are  necessary  but  the  seventh  one  simplifies  the  algorithm 
significantly.  The  method  was  based  on  the  solution  of  the  corresponding 
direct  problem  [2] .  Furthermore,  the  problem  of  acoustic  scattering  by  a 
soft  ellipsoid  coated  by  a  penetrable  confocal  ellipsoidal  shell  has  also  been 
solved  in  the  low-frequency  realm  in  [3].  The  present  paper  aims  towards  the 
development  of  an  algorith  that  will  provide  the  size,  the  orientation  and  the 
shell  thickness  of  a  soft  scatterer  which  is  coated  by  a  penetrable  shell  with 
exterior  boundary  an  ellipsoid  possessing  the  same  foci  as  the  scatterer 
itself.  It  is  shown  that  knowledge  of  the  first  two  low-frequency 
approximations  of  the  scattering  cross-section  for  specific  directions  of 
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iacidoce  ane  necesssay  to  identify  t&e  scattercr  ocofiletely.  Mare 
acesratoiy:  ve  need  to  ksew  tie  ^'a11Ee  of  tiie  leadis$  tom  for  asy  directicKs  of 
oxcitatiiKB  asd  t&e  valce  of  tbe  seco&d  tom  for  six  specific  directicMus  of 
excitation.  For  scalar  scattcria'  ly  aa  ellipsoid  o&e  can  look  is  [6,7,SJ . 
Inverse  scattcrlss  ty  an  ellipsoid,  from  a.  conpletely  different  point  of  viev, 
is  developed  in  [1]. 


1.2  _  ,2  2  _  ,2  ,2 
"3  "  ■  ®2  "  "  ^2 

Suppressing  the  harmonic  time  dependence  and  assuming  the  incident  plane 
wave 

i(r)  =  c'i-J,  (4) 

the  direct  scattering  problem  at  hand  asks  for  the  evaluation  of  a  function 
t  (r)  defined  on  the  domain  T  between  Sj^  and  and'  a  function  t'‘'(r)  defined 
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oa  liic  exterior  to  S_  dosaio  vhicli  satisfy  the  following  conditions.  The 

3 

total  field 


♦*(£)=♦(£)+%)»  (5) 

where  U(r)  is  the  scattered  field,  has  to  satisfy  the  Heliholtz  equation  with 
wave  snaiber  k 

(4«kV(jC)  =  0.  (6) 

while  the  inner  field  f (r)  has  to  satisfy  the  Helaholtz  equation  with  a 
different  wave  nuaber  K 

(4.kVe)  =  0,  rer.  (7) 

On  the  outer  boundary  the  boundary  conditions 
'a 

♦■(£)  =  Hi).  I  e  (8) 

(9) 

that  secure  penetration  are  to  be  satisfied,  while  the  inner  boundary  Sj^ 
consists  of  a  soft  boundary 


*-(r)  =0,  r€S|,  (10) 

The  positive  constant 


P 

B  =  — 

P' 

determines  the  ratio  of  the  outer  to  the  inner  mass  density, 
the  bulk  moduli  of  elasticity  is  given  by 


(11) 

The  ratio  of 


(12) 


where  ij  stands  for  the  relative  index  of  refraction.  For  convenience  we 
introduce  the  parameters 


B  =  B  -  1 
G  =  Bg^  -  1 
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(13) 

(14) 


which  imply 


,  C+1 

’  "  B+T* 


(15) 


Finally,  the  scattered  field  U(r)  should  satisfy  the  Somraerfeld  radiation 


condition  at  infinity. 

The  low-frequency  analysis  of  this  scattering  problem  can  be  found  in  [3]. 
In  particular  the  leading  two  low-frequency  coefficients  of  the  scattering 


cross-section  are  provided  by 


+  O(k^) ,  k-.0 


(16) 
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ijw 


1  r+00  dx 

"  2  J/.4  /— 2  r-2  / 
1  yxtaj  yx+a2  y: 

H"  =  B  iJ(aj)  +  iJ(bj) 


X+an 


(22) 


(23) 


and 


T  = 


o  L  loW'i'’  *  -;;;o  j, 


C(B+1) 
2t(H®)^  JV 


(afb2){Bri) 


5^[(-12(B+l)2t3(tol)(B+4)-B  ^].lj(a,)- (Bt2)^  iJ(bj) 


6(H") 


B+1t^  2C(B+1)  . 

^  (aj_a2a3-b^b2b3)  Ijj(bj) 


3(r)^ 


(24) 


2C(B+l)^aja2a3 


J  1x2 

“  D 


3(11°)^ 


J(ai)-lJ(bi))lJ(ai)-  ^[(2B.l)(B.l)-  iB(Cri)] 


2(b^+b2+b3)(B+l) 


3(11®)^ 


'{B2t2B-C)lJ(aj)*(C+l)lJ(bj)]-[— 


Btli^  bj+bj+bl 


^^^[(C+1)  E  bJl?(b,)n(B^+2B-C)  ^ 

3(H‘>)3|.  n=l  "  1'  1'  '  ■'„=1  ■' 


The  expression  for  T  involves  only  the  two  parameters  B,  C  and  the  six 

fV 

semiaxes  a^  a2,  ag,  b^,  b2,  b3.  Furthermore,  the  dyadic  B  corresponds  to  an 
ellipsoid  which  is  reciprocal  to  (2). 
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3.  THE  INVERSE  PROBIEH 


The  inverse  problem  we  consider  here  asks  for  the  determination  of  the 
semiaxes  a^  a2,  a^,  bj,  b2,  b^  which  specify  the  size  and  the  thickness,  and 
the  three  Euler  angles  B,  y  that  provide  the  orientation  of  the  scatterer. 
The  form  of  (16)  indicates  that  it  is  possible  to  use  the  method  developed  in 
[4]  to  solve  the  present  problem  as  well.  The  only  difference  is  due  to  the 
ellipsoidal  shell  that  coats  the  soft  scatterer,  i.e.  we  need  to  determine  bp 
b2,  bg  and  d,  y  as  in  [4],  as  well  as  the  semiaxes  a^  a2,  a^. 

Nevertheless,  because  of  the  confocality  of  and  Sj^,  relations  (3)  can  be 
used  to  determine  ap  a2,  in  terms  of  a^  as 

a?  =  ,  a|  -  b|  (25) 

,2  _  .  2  2  ,2  /0R\ 

3-2  *“  ^2  3^  *• 

Consequently,  we  only  need  to  determine  bp  b2,  bg,  ag  and  ji,  0,  y. 

An  expression  that  related  ag  with  bp  b2,  bg  can  be  obtained  by  a  single 
measurement  m^  of  the  leading  approximation  of  a. 

Indeed,  if 


rB+1-, 


.0  =  4i 


ll°J 


then,  in  view  of  (18),  (19)  and  (23)  we  obtain 
f+oo  dx 


B 


0 


^x+(bj-b|)+a|  yx+(b|b^)+a| 


(27) 


(28) 
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which  furnishes  as  a  function  of  b^,  b2,  bg 

a^  —  a0(bj)  ^2’  ^3^* 

The  existence  of  the  function  (29)  is  secured  from  the  implicit  function 
theorem  with  the  help  of  the  relation  [2] 


=  1 


(30) 


which  connects  the  corresponding  elliptic  integrals. 

Following  the  procedure  described  in  [4]  we  excite  the  scatterer  from  the 
six  directions 


k 


1’ 


k'  k' 


O 


xi+xl 

^  A// 


xA+Xo 


> 


(31) 


of  an  arbitrary  coordinate  system  indicated  with  primed  variables.  If  we 
denote  by  Xp  X2,  Xg  the  coordinate  system  that  fits  the  principal  axes  of  the 

/V 

target  ellipsoid  then  there  exists  an  orthogonal  dyadic  P  such  that 

r  =  P.r'  (32) 

A/ 

where  the  components  of  P  are  known  trigonometric  functions  [4]  of  the  three 
unknown  Euler  angles  (j>,  0,  y. 

Let 
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lrB+1 


2 


n.  =  4x  T  -  o — 7  k-.B.k. 

1  [  o[  j|OJ  ~i  ~i 


(33) 


be  the  measured  values  of  the  second  approximation  of  the  scattering 
cross-section  whenever  the  direction  of  excitation  is  given  by  k^ 
i  =  1,2, 3, 4, 5, 6.  A  series  of  calculations  similar  to  those  in  [4]  leads  to 
the  dyadic  equation 


<N/  12tT 

B  -  I  =  P.M.P* 


where  I  is  the  idemf actor,  P*  is  the  adjoint  to  P  dyadic  and 

3rA  ^  A  ^  AM 


(34) 


M= 


2m, 


■2niix^  2m2X2 

+  (mj+m2- 2m^) (xj  0  X2+X2  ®  x^) 

A  A  A  A 

+  (m2+ni3-2!n5)(x2®X3+X3(g>X2) 

A  A  A 

(m3+m^-2mg)(x3(x)xj+x^(x)x3)^ 
is  the  dyadic  of  measurements. 


(35) 


Since  M  is  a  real  symmetric  dyadic  its  normal  form  involves  three  real 
eigenvalues  and  an  orthogonal  set  of  eigenvectors. 

r-j  t\j 

By  virtue  of  (34),  the  eigenvectors  of  M  form  the  dyadic  P  while  its 
eigenvalues  Ap  A2,  A3  are  given  by 


K  =  K- 

n  n 


12xT 


n  =  1,.\3. 


(36) 


Hence  the  square  of  the  three  semiaxes  of  the  core  of  the  scatterer  are 
provided  by 


12xT 


b^  =  A  +  — 
n  n  m„ 


-,  n  =  1,2,3. 


(37) 


In  view  of  (29)  T  is  a  function  of  bj,  b2  and  b3.  Therefore,  (37)  forms  a 
highly  non-linear  system  of  three  equations  for  the  three  unknowns  b^  b2,  b3. 
Once  a  numerical  scheme  furnishes  the  value  of  b^,  b2  and  b3.  Formula  (29) 
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provides  Eg  and  (25),  (26)  give  the  values  of  a^^  and  Eg.  This  procedure 
identifies  the  size  of  the  scatterer  and  the  thickness  of  its  core.  The 
Euler  angles  that  give  the  orientation  as  well  are  obtained  from  the  knowledge 

A/ 

of  the  dyadic  P  exactly  the  same  way  as  in  [4] . 

We  mention  here  that  if  the  size  and  the  orientation  of  the  exterior 
ellipsoid  are  known,  then  a  single  measurement  of  the  leading  term  of  the 
scattering  cross-section  sufficies  to  identify  the  size  of  the  core  ellipsoid 
[5] .  This  procedure,  which  forms  a  nondistructive  method  of  evaluating  the 
size  of  the  interior  ellipsoid,  is  based  on  the  Rayleigh  approximation  of  a 
soft  ellipsoid 


and  the  corresponding  approximation 

.  ~  ]  (39) 

given  by  (16) . 

Eliminating  iQ^(aj)  between  (38)  and  (39)  we  obtain  the  relation 

4  rB+l  B  ■ 

I„(b()  =  Ur[jr  -  ,  (40) 

which,  for  known  a  and  c',  can  be  solved  numerically  to  obtain  bg  and  then  b^ 
and  b2  from 

bj  =  +  b|  (41) 

b^  =  hj  +  bj  (42) 

as  long  as  h2  and  hg  are  known  from  the  knowledge  of  Ep  a2,  Eg. 
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G.  DASSIOS,  K,  KIRIAKI  AND  V.  KOSTOPOULOS 

Inverse  thermoelastic  Rayleigh  scattering 
by  a  rigid  ellipsoid 

ABSTRACT 

In  this  paper  the  inverse  scattering  problem  for  the  rigid  ellipsoid  in  linear 
thermoelasticity  is  examined.  We  prove  that  six  measurements  of  the  far 
field  pattern  in  the  low-frequency  region  are  necessary  in  order  to  evaluate 
the  semiaxes  of  the  ellipsoid  as  well  as  to  fix  the  position  of  its  principal 
axes. 

1.  INTRODUCTION 

The  inverse  scattering  problem,  as  it  is  well  known,  is  concerned  with  the 
problem  of  determining  the  shape  and/or  the  physical  properties  of  the 
scattering  object  from  the  knowledge  of  the  scattered  far- field  data. 

The  mathematical  methods  used  to  investigate  the  inverse,  as  well  as  the 
direct,  scattering  problems  depend  heavily  on  the  frequency  of  the  incident 
wave  [11,12]. 

In  the  low-frequency  region  there  is  the  problem  that  low  frequency  data 
does  not  provide  enough  information  for  a  sharp  resolution  of  the  scattering 
surface  and  the  optimisation  procedure  requires  a  direct  scattering  problem  to 
be  solved  at  each  step  of  the  iterative  scheme  for  arriving  at  a  solution. 

In  low  frequencies  Angell  and  Kleinraan  [1]  described  a  method  for  finding  the 
dimensions  and  orientation  of  an  ellipsoid  from  a  constrained  optimisation 
problem  for  a  functional  defined  in  terms  of  the  polarisability  tensor 
elements  associated  with  the  object.  Dassios  [5]  solved  the  inverse 
scattering  problem  for  the  acoustical  soft  ellipsoid  based  on  the  informations 
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which  the  direct  problem  provides  in  the  radiation  region.  The  inverse 
scattering  problem  in  linear  elasticity  was  examined  in  [2] . 

In  addressing  oneself  to  the  practical  problem  of  reconstructing  the  shape 
of  an  obstacle  from  far- field  data,  one  is  faced  with  problems  of  numerical 
instability.  The  inverse  scattering  problem  in  general  is  ill-posed,  in  the 
sense  of  Hadamard,  and  intrinsically  nonlinear.  So,  the  aim  of  the  provided 
methods  for  solving  inverse  scattering  problems  is  to  be  relatively  easy  in 
numerical  implementation.  A  survey  of  the  research  done  in  this  area  can  be 
found  in  [3,4,13]. 

The  inverse  scattering  problem  which  is  considered  in  this  paper  is  to 
determine  the  semiaxes  and  to  fix  the  orientation  of  a  rigid  ellipsoidal 
scatterer  embedded  in  an  infinite,  homogeneous,  isotropic  thermoelastic  medium 
from  the  knowledge  of  the  far- field  data,  in  a  finite  number  of  directions. 

In  Section  2  we  formulate  the  direct  scattering  problem  and  we  give  the 
necessary  results  obtained  in  [7,8] . 

In  Section  3  we  examine  the  inverse  scattering  problem  for  a  rigid 
ellipsoid  in  linear  thermoelasticity.  We  show  that  the  necessary 
measurements  in  order  to  obtain  all  the  information  needed  about  the 
ellipsoidal  scatterer  are  six.  This  is  a  consequence  of  the  simplicity  of 
the  expressions  for  the  leading  terra  of  the  angular  scattering  amplitudes  in 
low-frequency  regions  that  allows  for  an  exact  solution  of  the  corresponding 
inverse  problem.  The  nonlinearity  of  the  problem,  expected  for  inverse 
problems,  enters  via  the  elliptic  integrals,  the  second  power  of  the  values  of 
the  semiaxes  and  the  quadratic  expressions  of  the  components  of  the  directions 
of  incidence. 

Finally,  in  Section  4  remarks  for  the  stability  of  the  numerical  method 
and  a  general  discussion  concerning  our  approach  are  presented. 
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2.  FORMULATION  OF  THE  DIRECT  SCATTERING  PROBLEM 


Let  us  assume  that  the  solid  ellipsoid 

3 

S  “2  <  Ij  0<a2<a2<aj<+oo  (1) 

i=l  a^ 

is  embedded  in  an  infinite,  homogeneous  and  isotropic  tkermoelastic  medium. 

The  isotropic  thermoelastic  material  is  characterised  by  the  Lame 
constants  A,/i,  the  mass  density  p  the  coefficient  of  thermal  diffusivity  a, 
and  the  linear  expansion  coefficient  a.  The  unified  four- dimensional  field 

U(l)  =  (ui(r),  U3(r),  B(r)),  (2) 

where  u(r)  denotes  the  displacement  and  0(r)  the  temperature  field,  specifies 
the  stationary  thermoelastic  state  of  the  medium  whenever  it  belongs  to  the 
kernel  of  the  time- independent  Biot  operator 


a; 

and  Ijj  denotes  the  unit  dyadic  in  n- dimensions.  A  convenient  dimensionless 
coupling  constant  is  provided  by  the  parameter 

^  ""  'MJi 

In  the  limit  as  7  ->  0+,  ?/  -»  0+,  e  -»  0+  the  thermoelastic  problem  decouples  to 
the  corresponding  scattering  problem  in  the  classical  theory  of  elasticity  and 
an  independent  heat  conduction  problem. 

A 

A  general  incident  plane  wave  propagating  in  the  direction  k  has  the  form 

[7] 
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where 


(6) 


t(i)  =  +  i^(i)  +  t*’(r) 


(7) 


(8) 


(9) 


1  ^  -  ikA.r 

r(r)  =  Al(k,/?j)e  - 

is  the  elastothermal  plane  wave  of  amplitude  A^. 

n  A  ikrtk  •  r 

r(r)  =  A^(ik2/^2J^,ik2)e  ‘ 

2 

is  the  thermoelastic  plane  wave  of  amplitude  A  ,  and 

^  -  ik^k.r 

l®(r)  =  A^(b,0)e  ■ 

is  the  transverse  plane  wave  of  amplitude  A®  and  polarisation  along  the 

direction  b,  orthogonal  to  the  direction  of  propagation  k.  The  factor  ik2  in 

(8)  has  been  added  in  order  to  secure  analyticity  of  §  with  respect  to  the 

Is  9 

wave  number.  The  amplitudes  A  and  A  have  dimensions  of  length,  while  A 

has  dimensions  of  length  times  temperature. 

In  consistency  with  physical  reality,  the  wavenumbers  kp  k2  are  chosen  to 

be  those  roots  of  the  characteristic  systems  [7] 


^2  ■  9(l+e)  +  kp 


k^k^ 


=  qk: 


(10) 


for  which  Imk.  >0,  i  =  1,2.  These  conditions  reflect  the  dissipative 
character  of  the  thermoelastic  medium.  The  constants 

kjW^ 

^1  "  ■  T2~ ’ 


kj-q 
“  (k2-k2)(A+2„) 


(12) 


1  9 

furnish  the  appropriate  factor  in  order  for  and  to  belong  to  the  kernel 
of  L.  On  the  other  hand,  as  e  0+,  -»  0  and  /^2  9  and  the  system 
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(13) 


^ecoEplcs  Ib  a  asteral  way. 

Ta  classical  clasticiiyx  w'C  liavc  tlsc  wave  relacimis 

fir  =  k  c,  =  c 
p  p  s  s 

wkere  aad  k^  dcaotc  tiic  wai'caieiiicrs  of  tkc  clastic  loisitndiaal  and 
traasvcrsc  waves,  aad 


^2  ^  A*2p  2  f* 

P  "  ~P~* 

specify  the  correspoadiB'  pkasc  \'clocitics. 

In  thciaoclasticity,  vc  denote  t&ic  coo^lcx  wavemmaers  kj  and  k.^  by 


a 

1  =  57  ^  “p 

Vj  >  0, 

■>1 

>  0 

V2  >  0, 

^2 

>  0 

(14) 


(15) 


where  v^.  V2  are  the  phase  velocities  of  the  elastothcmal  and  themoelastic 
waves,  respectively,  and  dp  d2  deteroine  the  corresponding  dissipation 


coefficient.  In  the  decoupled  case,  as 

£  -!  O-j- 

kl-‘p  1  (>6) 

kj  ->  /4  =  (W)  7s  i 

The  general  scattered  field  U  has  a  corresponding  dccoaposition  into 
clastothennal,  thcrmoelastic  and  transverse  parts  via  the  relations 

E(r)  =  n'(r)  +  a^(r)  +  n®(r),  jj-j 

a(r)  =  e‘(r)  t  B^Cr)  ) 

A  set  of  ten  asymptotic  relations  given  by  Kupradze  express  the  radiation 
conditions  which  secure  the  well-poscdness  of  our  scattering  problem. 

The  unified  total  field  is  expressed  through 

i(r)  =  i(r)  +  U(r)  (18) 

The  boundary  conditions  for  the  rigid  ellipsoid  in  thermoelasticity  are 
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descri!»cd 


(19) 


«l(ar.l!)*(d  =  4.  E  e  S.  k  =  >=2. 


vherc  the  bonndary  diffcreatial  operators  are  expressed  as  follows: 
(i)  For  the  rigid  ellipsoid  at  zero  teaperature 


h 

0 

0 

1 1 J 

B,(aj.,n)  =14  = 

(ii)  For  the  theraa^ly  insulated  rigid  ellipsoid 


(20) 


I3 

0  ■ 

0 

(21) 


n  being  the  outward  unit  nomal.  The  integral  postulation  of  the  above 
problca  is  presented  in  [7] .  The  far- field  behaviour  of  the  scattered 
elastic  and  theraal  fields  have  been  studied  in  [7]  where  the  following  six 
noraalised  scattering  auplitudes  have  been  introduced. 


As  r  -1  4»,  the  following  asymptotic  forms  hold  true 

-d.r 


u\r)  =  e  ^  [gj(r,k)h(—  r)r  +  0(r'^) 


(22) 


ll^(r)  =  e  ^  [gj(r,k)h(—  r)r  +  0(r  ^)] 


(23) 


eVl)  =  rji  +  OCr-2)] 


(25) 
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8^(r)  =  e  *  J<*(r,k)h(—  r)r+  0(r  ^)J 


(26) 


where  ■  =  1,2  define  the  radial  elastothennal  (for  m  =  1)  and  the  radial 
theraoelastic  (for  ■  =  2)  nomalised  scattering  aaplitude,  g°,  sl  =  0,  ^  define 
the  angular  normalised  scattering  amplitudes,  and  f”,  m  =  1,2  define  the 
corresponding  thermal  normalised  scattering  amplitudes  for  the  elastothermal 
and  thermoelastic  waves  respectively.  The  analytic  expressions  of  the 
thermoelastic  scattering  amplitudes  are  given  in  [7] . 

In  the  theory  of  thermoelasticity,  where  the  existence  of  dissipation 
manifests  itself  via  the  nonvanishing  imaginary  parts  of  the  wavenumbers 
and  k2,  the  unit  of  energy,  which  enters  the  definition  of  the  scattering 
cross-section,  is  not  any  longer  invariant  under  translations  in  the  direction 
of  propagation. 

This  is  the  reason  for  introducing  in  [7]  a  local  determination  of  the 
scattering  cross-section  which  in  general  has  the  form 

e,e(k)  e“(k)«22(k).e|={k).el2jk) 


<r  = 


'in 


22^  ss^^l2 
e-  +e-  +e-  +6- 
in  in  in  in 


(27) 


where  e|^ 


-22 


in 


3?" 

in 


in 


is  the  part  of  the  unit  of  energy  that  corresponds  to  the 
elastothermal  incident  wave, 

is  the  part  of  the  unit  energy  that  corresponds  to  the  thermoelastic 
incident  wave, 

is  the  part  of  the  unit  energy  that  corresponds  to  the  transverse 
incident  wave, 

is  the  part  of  the  unit  of  energy  that  corresponds  to  the 
interaction  between  the  elastothermal  and  thermoelastic  incident 
waves , 
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and  ejl,  ell,  ell,  ell  describe  the  corresponding  parts  of  the  energy  of  the 
SC  so  so  so 

scattered  wave. 

From  the  analytical  expressions  of  the  above  quantities,  given  in  [7] ,  it 
is  easy  to  observe  that  only  the  terms  e^®,  e|®  are  independent  of  r,  while 
all  the  other  quantities  die  out  as  r  -»  +».  Therefore,  the  quantities 
®sc’’®in’  ®sc’  ®in’  ®sc  ®in  leading  role  in  the  definition  of  the 

scattering  cross-section,  which,  as  r  -»  +»,  has  the  asymptotic  form 
e®®(k) 

^in  s'  '  ‘ 

In  other  words,  the  scattering  cross-section  is  measured,  to  a  first 
approximation,  by  the  only  nondissipative  part  of  the  incident  and  scattered 
wave  which  survives  far  away  from  the  scattering  region,  i.e.  the  transverse 
part. 

CO 

For  the  above  equivalent  thermoelastic  scattering  cross-section  (r  ~  tr  , 
an  angular  type  scattering  theorem  as  in  classical  elasticity  [10]  is 
obtained,  from  which  we  have  the  scattering  relation 

In  low-frequency  region  the  leading  terms  of  the  angular  thernioelastic 
scattering  amplitudes  for  the  two  problems  that  will  be  examined  in  this  paper 
are  given  by  [8] 

g®  =  -  4^  a.J  [t' 11^(1' )Jds(r')  +  O(k^),  k  -»  +0,  a  =  0,(l>  (30) 

s 

where  for  the  case  of  a  rigid  ellipsoid  u^  is  the  solution  of  the  following 
boundary  value  problem 

T-pAUoCr)  +  (l-rp)VV.Up(r)  =0,  p  > 

]lo(r)  =0,  p  = 
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(31) 

(32) 


u^(r)  =  (A^  +  k\)i  +  k%  +  Oy ,  p  -.  +«,  (33) 

which  can  be  solved  by  introducing  a  fictitious  potential  in  the  Papkovich 
representation  and  then  use  the  freedom  it  provides  to  decompose  the  general 
boundary  condition  into  two  simpler  ones.  Following  the  method  introduced  in 
[6]  we  obtain,  after  very  long  calculations,  the  solution 

i„(r)  =  [(a‘  +  k\)k  +  A%]  . 


3 

.  S 

m=l 


/  2  2\ 
(/>  -«i) 


/  2  2  /  2  2 
y/  p  yj  p-V 


Q  (2l  x„  r 

o  -  m  ^  -  ra  ~ 

where 


(34) 


du 


/  2  2  2  /  2  2  2 
V  u  -  V  u  -  0'i+ff3 


(35) 


'+00 

du 

p  f  2  2  2\ 

7222/222 
yu  -0^+02  y « - 

(36) 
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(37) 


=  Il  =  l!(“l).  "  =  1.2.3 

and  g  is  the  outward  unit  normal  on  the  ellipsoid. 

The  expressions  of  the  scattering  amplitudes  involve  the  integral 

[  T  u  (r)  ds(r) 

where 

T  UqU)  =  2/i^.VuQ(r)  +  A£V.UQ(r)  +  Mx(VxuQ(r))  (38) 

provides  the  traction  field  on  the  surface  of  the  ellipsoid.  In  view  of  the 
potential  integral 


1 


ds 


=  4t 


'p=a^  /  2  2  /  2  2 

^  1  V  P  -/^  ^  p-v 
and  a  series  of  vector  and  dyadic  calculations  we  obtain 


(39) 


3  - 


(40) 

' p=a^  ''  *  '‘'p  *''“m*l  ''■p' 

Introducing  (40)  into  Eq.  (30)  we  conclude  that  the  angular  amplitudes  assume 
the  form 

3  -- 

'p  ''ml'p  '0 


I  1  .  *273 


3  [(AVA‘?2)VA\] 

Ji  (Ai)A»(Ai)ii  *  “(' ) 

'p  'ml'p  '0 


(42) 


So  from  Eq,  (28)  by  substitution  of  the  angular  scattering  amplitudes  and 


supposing  only  S-wave  incidence  we  conclude  that 


a 


qo_  3  ^ 

,ss  S  '■ 


3  m=l 


♦  0(k2) 


(43) 
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(44) 


where  b  =  (b|,b2,b2),  b.k  =  0  and 

■■S  =  (“l)  "  =  1>2,3 


3.  THE  INVERSE  SCATTERING  PROBLEM 

We  assume  that  the  directions  of  the  principal  axes  of  the  ellipsoid  described 
by  Eq.  (1)  coincide  with  the  unit  vectors  of  the  orthogonal  cartesian 
system. 

We  choose  six  arbitrary  directions  of  incidence  of  transverse  waves. 

This  choice  of  S- waves  is  dictated  by  the  above  discussion  about  the 
dissipation  of  all  but  e?®  and  e|®  in  the  far  field  region. 

In  the  sequel  we  will  see  why  six  measurements  are  enough  to  specify  the 
three  semiaxes  as  well  as  the  three  Euler  angles  that  fix  the  position  of  the 
principal  axes  of  the  ellipsoid. 

For  S- incidence  along  the  k- direction  we  can  evaluate  the  angular 
thermoelastic  scattering  amplitudes,  so  that  Eqs.  (28,  29,  43)  provide  the 
thermoelastic  scattering  cross-section. 

Let  the  six  arbitrary  directions  of  incidence  be 


k' .  =  x' 


j  =  1,2,3 


ki  =  —  (xA+xI) 

-5  '-2  -3' 


(45) 


where  {x|,K2  ,213}  form  an  orthonormal  set  of  an  arbitrarily  chosen  cartesian 
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r 


system  whose  origin  coincides  with  the  centroid  of  the  ellipsoid.  In  the 
sequel  with  primes  we  will  refer  to  the  arbitrarily  chosen  system  and  without 
primes  to  the  system  determined  by  the  principal  axes  of  the  ellipsoid. 

/V 

Let  P  be  the  orthogonal  matrix  that  transforms  (by  rotation)  the  system 

AAA  AAA 

5X23^3  system  XpX2>^3  which  determines  the  principal  directions  of 

the  ellipsoid,  that  is  we  have  the  relation 

fsj 

r  =  Pr' .  (46) 

The  knowledge  of  the  matrix  P  will  provide  the  orientation  of  the  ellipsoid. 

r>j 

The  elements  P- j  of  P  are  expressed  via  the  three  Euler  angles  (i,  6,  y  by  the 
relations 

Pjl  =  cos  cos  y  -  cos  6  sin  ^  sin  y 
P^2  =  sin  cos  y  +  cos  9  cos  ^  sin  y 
Pj2  =  sin  9  sin  y 

P22  =  -cos  ^  sin  y  -  cos  ^  sin  ^  cos  y 

P22  =  -sin  ^  sin  y  +cos  9  cos  ^  cos  y  (47) 

P23  =  sin  9  cos  y 

Pgj  =  sin  9  sin  ^ 

P22  =  -sin  9  cos  ^ 

Pgg  =  cos  9 

Since  the  relation  (43)  is  referred  to  the  principal  axes  system  it  follows 

r-j 

that  (43)  holds  true  after  the  transformation  described  by  the  matrix  P,  has 
been  applied  to  the  directions  of  polarisation  b',  that  is 

bj=Pbl  (48) 

So,  the  six  measurements,  in  view  of  Eq.  (43),  yield 

32t  ..  m 

nij  =  ^  -j’  j  =  (49) 
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where 


A.j^  =  0,  k 


(50) 


.  _  ^  i  =  k  =  1,2,3 

^ik  “  rl 
^0 


From  Eqs.  (48,  49)  we  conclude  that 


■j  =  f 


(51) 


Substituting  Eq.  (45)  in  Eq.  (51)  we  obtain  the  following  system  of  equations 
for  the  unknown  quantities  (iS,^,y,apa2,a2 

3 

A  .n  ■  32i  ’  ^ 

n=l 


(52) 


3('’nAj)  H 


n= 


=16;,  (i.j.k)  =  {(1,2, 4), (2, 3, 5), (3, 1,6)) 


From  Eqs.  (52)  we  conclude  that 

3  Pnifnj  3(2l»|,-nj-.p 


n=l  L' 


32x 


(53) 


0 


where  (i,j,k)  as  in  Eq.  (52). 

So,  we  have  the  matrix  relation 

A/  /V  a; 

A  P  =  M 


(54) 


where  M  is  a  real  symmetric  matrix  with  known,  from  the  six  measurements, 
elements  given  by  the  relations 

3m 


^ii  ~  16t  ^ 


(55) 


Hi-  =  M.,  = 


3(2mj^-m.-mj) 


327r 
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(i,j,k)  as  in  Eq.  (52). 

From  Eq.  (54)  we  conclude  that 


t\f  A/  A/  fsjfj 

A  =  P  M 


(56) 


Since  P  is  an  orthogonal  matrix,  M  is  a  real  symmetric  one,  and  A,  from  Eq. 

«v 

(50),  is  a  diagonal  matrix  it  follows  that  the  eigenvalues  of  M  are  given  by 

^  =  4’  n  =  1,2,3  (57) 

^0 

~T 

while  the  columns  of  P  are  the  corresponding  orthonormal  eigenvectors.  So, 

/V 

the  eigenvectors  of  the  known  matrix  M  specify  the  orientation  of  the 
ellipsoid.  In  order  to  evaluate  the  semiaxes  of  the  ellipsoid  we  apply  the 
following  procedure. 

From  Eq.  (57)  we  have  that 


1 


^  “  1  >2,3 

°  ''n 

From  Eq.  (44)  and  the  well-known  formula 

i  ah’!  =  i' 

«=1  ”  1  0 

which  relates  the  elliptic  integrals  we  derive  that 

H  13  1 

ii(a.)  =  — 2—  ^  r  = 

®  ^  2(t“+2)  n=l  ''n 

From  Eqs.  (44,  57)  we  conclude  that 


(68) 


(59) 


.2Tn/ 


rhl 

1  1  V  3  1 


^n^l(^l)  "  — 2 -  ■  - 2 -  ^  “  1>2>3  (60) 

^  ^  (r^-1)  ''n  2(r^2)  i=l  ''i  « 

where  M^,  are  known  quantities. 


1 

In  order  to  bring  the  elliptic  integral  1^  to  its  canonical  form  we 
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perform  the  transformation 


o2  J 
a^-ag 


X  = 


(61) 


and  we  conclude  in  standard  notation  of  elliptic  integrals  that 

dt 


•  2 

t  sin  a„ 


(62) 


where 


(^0  = 


2  2 


(63) 


Oq  =  sin 


2  2 

"2^ 

“l-"3 


In  order  to  bring  the  elliptic  integral  of  the  second  kind  I^  to  its  canonical 
form  E(j5Q,aQ),  we  perform  the  same  as  above  transformation,  while  for  the 
integral  l|  we  conclude 


t^dt 


^2  .  2 
t  sin  a 


=  -(a2-a2,3/2lJ(a,)  (64) 


0 


Finally  after  some  algebraic  calculations  we  conclude  that 

ll(ai)  =  (af  a2)-3/2  [e((S  a  )-r(((  a 

Sin  a^  0''-' 

If  we  introduce  the  notation 


(65) 


■  a. 


■  aj 


(66) 


use  the  relations 
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1 


3 

S  I 
n=l 


n 

1 


\  h  -  -^0 

n=l 

and  the  equations  (62)  and  (64)  we  conclude  the  relations 


2  2„  2u  2u 

sp2M^+S2M2+SjMg 

®1®2 


(l-s2)-V2F(^^,a^)  =  M.  (67) 


SjS2M^+S2M2+s^M^ 


hh 


where  in  view  of  (66) 


(68) 


(l-s|)-l/^l-s2)-'[E(#„.a„)-F(^^,a„) 

=  Aj  (known  quantity) 

W  sin-'(l-s^)’/2 
1.w2 

,  ‘  1/2 

\  (73) 

I-S2 

In  order  to  solve  numerically  the  nonlinear  system  given  by  equations  (67,  68) 
we  use  the  same  iterative  scheme  as  in  [2] . 


(69) 


4.  DISCUSSION 

In  this  paper  the  inverse  scattering  problem  of  linear  thermoelasticity  for  a 
rigid  ellipsoid  in  the  low  frequency  region  is  examined.  Ve  use  information 
from  the  solution  of  the  direct  thermoelastic  scattering  problem  and  propose  a 
method  for  the  corresponding  inverse. 

In  order  to  evaluate  the  dimensional  characteristics  and  the  orientation 
of  the  scatterer  we  need  six  measurements  of  the  far- field  data  in  the 
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low-frequency  region  under  the  condition  that  we  have  a  knowledge  of  the  shape 
and  the  boundary  a  priori  conditions.  This  type  of  approach  of  the  inverse 
problem  has  been  proposed  in  [5] .  With  a  similar  technique  the  scattering 
problem  in  linear  elasticity  has  been  solved  in  [2] . 

If  the  boundary  conditions  on  the  surface  of  the  scatterer  alter  the  above 
approach  c?n  also  be  used.  Obviously  from  the  solution  of  the  direct  problem 
we  obtain  a  highly  nonlinear  system  which  cannot  be  solved  by  a  simple  and 
rapidly  convergent  iterative  scheme. 

Looking  closer  at  the  results  contained  in  this  work,  one  can  see  no  basic 
difference  if  we  consider  a  thermally  insulated  rigid  scatterer  or  a  rigid 
scatterer  at  zero  temperature.  In  other  words,  the  results  are  independent 
of  the  temperature  boundary  condition  on  the  surfac(  of  the  rigid  scatterer. 

This  is  a  reflection  of  the  fact  that  thermal  effects  do  not  enter  the 
Rayleigh  approximation  of  the  transverse  field  in  the  radiation  zone. 

The  corresponding  temperature  field  appeared  in  the  low-frequency 
approximations  of  order  higher  than  the  leading  one.  Even  in  these  higher 
order  approximations  the  dependence  on  the  temperature  is  implicit  through  its 
effects  on  the  particular  form  of  the  displacement  field. 
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A.FARIDANI 

Reconstructing  from  efficiently  sampled 
data  in  parallel-beam  computed 
tomography 

1.  IKIlflDUCriHf 

In  CG3ipQtcd  toanograpliy  (CT)  an  object  is  exposed  to  radiation  vhicb  is 
Bccasared  after  passing  through  this  object.  Fnn  these  ncasareisents  a 
certain  fanction  f ,  characterizing  the  interaction  of  the  object  with  the 
radiation  is  rcconstnictc-d.  This  function  represents  an  'iaage*  of  the 
interior  of  the  object. 

In  diagnostic  radiology,  the  classical  application  of  CT,  X-rays  arc  used 
as  radiation.  The  ■casurenents  are  then  line  integrals  of  the  X-ray 
absorption  coefficient  f .  This  leads  to  the  ■athexatical  problea  of 
reconstructing  a  function  froa  its  line  integrals.  In  tvo  disensions  this 
eeans  the  inversion  of  the  Radon  transfora. 

In  this  paper  ve  exaaine  how  nany  line  integrals  have  to  be  aeasured  in 
order  to  achieve  a  certain  accuracy  and  resolution.  Ve  confine  our 
investigation  to  the  so-called  parallel- beam  sampling  geometry.  First  we 
describe  the  application  of  multidimensional  sampling  theory  to  sampling  the 
two-dimensional  Radon  transform.  This  approach  was  first  taken  by  Lindgren 
and  Rattey  [11,  16]  and  further  developed  by  Natterer  [12].  It  leads  to 
sampling  schemes  which  need  a  minimal  amount  of  data  to  ensure  that  the  Radon 
transform  Rf  is  determined  by  the  measured  values  up  to  a  small  error. 

The  question  of  how  to  achieve  good  reconstructions  from  such  data  has 
been  studied  by  Kruse  [9]  who  obtained  error  estimates  for  the  filtered 
backprojection  algorithm,  the  most  popular  reconstruction  method.  The  main 
purpose  of  this  paper  is  to  derive  estimates  extending  Kruse's  results  and  to 
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use  tliis  theory  to  obtain  iaproved  rcconstmctions.  The  paper  is  organized 
as  follovs: 

In  the  remainder  of  this  section  ve  give  a  brief  description  of  the  Radon 
transfom  and  its  inversion.  In  the  next  section  ve  derive  a  version  of 
Petersen  and  Middleton's  saopling  theoren  [i4]  suitable  for  sanpling  functions 

defined  on  [0i2r]  x  t.  .  Our  proof  differs  fron  the  approach  taken  in  [12. 
p.64]  and  is  based  on  a  Poisson  suaaation  foniula  derived  in  [9] . 

Furthemore  ve  characterize  the  class  of  sanpling  schcncs  to  which  the 
sanpling  theorca  is  applicable  and  which  arc  suitable  for  sanpling  the 
two-dinensional  Radon  transfom  in  practical  applications.  Vc  call  these 
schcacs  'adnissiblc  sanpling  lattices*.  Application  of  the  sanpling  theorca 
leads  to  sanpling  conditions  which  have  to  be  satisfied  by  a  sanpling  schcoc. 
The  results  of  [11,  12,  16]  for  the  usually  enployed  standard  lattice  and  the 
so-called  'interlaced  lattice*  arc  derived.  The  interlaced  lattice,  first 
suggested  for  sanpling  the  Radon  transfom  by  Connack  [2] ,  requires  a  minimal 
ancunt  of  data.  Sampling  conditions  for  the  standard  lattice  were  already 
derived  in  different  ways  by  Bracewell  [1],  and  by  Crowther,  De  Rosier  and 
Klug  [3]. 

To  clarify  the  question  of  accuracy  of  reconstructions,  we  give  a  detailed 
error  analysis  for  a  particular  reconstruction  method,  the  filtered 
backprojection  algorithm.  In  section  3  we  describe  the  implementation  of  the 
algorithm  for  data  sampled  on  an  admissible  sampling  lattice.  In  the 
following  section  we  derive  estimates  for  the  reconstruction  error.  Ve 
extend  the  results  of  [9]  by  taking  into  account  the  influence  of  a  certain 
interpolation  step  occuring  in  the  algorithm.  This  interpolation  was 
neglected  in  [9].  It  turns  out  that  it  is  critical  when  the  interlaced 
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lattice  is  used,  but  that  good  reconstructions  can  be  achieved  with  a  sensible 
choice  of  the  paraaeters  of  the  algoritha. 

The  proof  of  the  estimates  already  derived  by  Kruse  is  closely  related  to 
the  proof  of  the  saapling  theorea  given  in  section  2.  This  indicates  that 
the  filtered  backprojection  algoritba  aight  be  aore  suitable  than  other 
reconstruction  aethods  when  ainiaal  saapling  scheaes  like  the  interlaced 
lattice  are  used. 

The  last  section  is  devoted  to  nuaerical  experiaents  and  their  discussion 

in  the  light  of  the  results  of  section  4.  It  turns  out  that  the  theoretical 

results  explain  the  iaainent  nuaerical  difficulties  and  show  the  way  to  obtain 

good  reconstructions,  thereby  improving  previous  numerical  results . 

In  the  following  we  give  a  brief  introduction  to  the  Radon  transform  and 

its  inversion  and  also  introduce  some  notation. 

Let  R,  N,  Z  denote  the  real  numbers,  natural  numbers,  and  integers, 

respectively.  Furthermore  let  Q  denote  the  unit  disk  in  R  and  T  the 

interval  [0,2x).  For  (5  e  T  the  variable  ff  will  always  denote  the  unit  vector 
T 

(cos  (/,sin  (})  . 

The  Radon  transform  of  a  function  f  e  C^(fi)  is  given  by 

Rf(^,s)  =  r  If  (s^+t^dt,  6  T,  s  e  R.  (1) 

J-00 

For  a  survey  on  the  mathematical  properties  of  the  Radon  transform  and  its 
many  applications  see  e.g.  [4,  6,  7,  8,  12]. 

As  we  will  see  in  Theorem  1.1  below,  the  Radon  transform  is  closely 
related  to  the  Fourier  transform.  The  Fourier  transform  of  a  function  F  6 
Lj(R’')  is  given  by 

KO  =  [  F(x)e"^'^-^dx  (2) 

where  x.^  =  x^^^  denotes  the  dot  product.  The  Fourier  transform  can  be 
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extended  to  the  space  of  teapered  distributions.  Its  inverse  is  given  by  the 
fonnla 

F(x)  =  F(-x).  (3) 

The  convolution  f*g  of  two  L.-functions  is  given  by  f*g(x)  =  f (y)g(x-y)dy. 

Ve  have 

=  (2x)”/^i(Og(f)-  (4) 

In  case  of  the  two-dimensional  Radon  transform  the  Fourier  transform  with 
respect  to  the  second  variable  alone  is  also  of  interest.  It  is  defined  by 

f  M((i,s)c-*'^ds  (5) 

”2 

where  indicates  that  the  Fourier  transform  is  taken  with  respect  to  the 
second  variable  only. 

The  Radon  transform  is  closely  related  to  the  Fourier  transform  by  means 
of  the  following  theorem: 

Theorem  1.1 

LetfGC^(fl).  Then 

(Rf)'2{d,»)  (6) 

Proof:  Replace  Rf((},s)  in  (5)  by  the  righthand- side  of  (1)  and  use  the 
definition  (2)  of  the  Fourier  transform.  □ 

Taking  a  two-dimensional  inverse  Fourier  transform  on  both  sides  of  (6) 
yields  immediately  an  inversion  formula  for  Rf : 


=  l(2ry^  r 

JQ  -*-00 

=  5(2z)-*/2  f'  r  l,|(M)'2((l,,r)e'“  ‘'drf!i 
^  Jo  J-OO 

(7) 

1  f2x 

=  4rJ^ 

d 

where  q((5,s)  =  Rf((J,s),  and  H  denotes  the  Hilbert  transform  acting  on  the 

second  variable.  For  other  types  of  inversion  formulas  see  e.g.  [12,  Chapter 
II]. 

It  is  readily  seen  from  (7)  that  an  exact  inversion  of  the  Radon  transform 
is  unstable  because  of  the  amplification  of  high  frequencies  due  to  the  filter 
\a\  in  the  inner  integral  of  (7).  More  suitable  for  numerical  inversion  are 
approximate  inversion  formulas  of  the  following  kind,  where  this  instability 
is  removed  by  means  of  a  suitable  low- pass  filter: 

Theorem  1.2 

2 

Let  WjjCR  -»  R  be  a  radially  symmetric  low- pass  filter  with  cut-off 
frequency  b,  i.e. 

\(.0  =  {2')‘'«lfl/b) 

with  0  <  <  1  and  ^((7)  =  0  for  \a\  >  1.  Define  Wj^  by 

and  let  f  e  C^) .  Then 

2t 

=  J  J  Wjj{x.0'S)Rf (j5,s)dsd<}.  (8) 

Proof:  Proceeding  as  above  we  obtain 

W|,*f(x)  =  [ 
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1  At  rflo 

L 


kll5(kl/b)f((r^)e^‘^-^d<rd(5 


=  j  j  W|j(x.^s)Rf((},s)dsdj} 


where  we  have  used  (4).  d 

This  kind  of  approximate  inversion  formulas  goes  back  to  [19] ,  and 
provides  the  basis  for  the  filtered  backprojection  reconstruction  algorithm 
[15] ,  which  we  will  describe  in  section  3.  There  are  many  other  methods  for 
the  numerical  inversion  of  the  Radon  transform.  For  a  survey  of  such 
reconstruction  methods  see  [12,  Chapter  V]. 

In  tomography  one  has  to  compute  approximations  to  f  from  measurements  of 
Rf .  The  question  arises,  how  many  measurements  are  needed  and  at  which 
points  ({},s)  we  should  sample  Rf({},s).  Therefore  we  now  discuss  the  relevant 
sampling  theorem. 


2.  SAMPLING  ON  [0,2t)"^  x 

Since  Rf({},s)  is  2t- periodic  in  the  first  variable,  we  need  a  sampling  theorem 
for  functions  which  are  2t- periodic  in  some  of  their  variables.  They  can  be 

^1  ^2 

regarded  as  functions  on  T  x  R  ,  where  T  denotes  the  interval  [0,2t).  One 
essential  tool  for  deriving  the  sampling  theorem  is  the  Fourier  transform  on 

T  X  R  which  is  defined  as  follows:  Let  S  ^  ^  denote  the  class  of 
(f- functions  of  n^  +  n2  variables  which  are  2t- periodic  in  their  first  n^ 


iirt 

variables  and  in  C^(R  )  with  respect  to  the  last  n2  variables.  The  Fourier 


,n 


1  ^^'9 

transform  of  a  function  G  G  S  is  a  function  G:Z^xR^-iC  given  by 


1 
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=  f,  CW,s)e-»''-V5k-<'dsd^ 

T  ^  R  ^ 

n^i  Ha  n*!  Hrt 

with  (k,<r)  G  Z  ^  X  R  ^  and  (^,s)  €  T  ^  x  R 

Let  Lj(Z  ^  X  R^)  denote  the  space  of  all  functions  g  =  Z  xR^-dC  for 


which 


^  [  n  <  «>• 

tl  .  J  ^ 


kGZ 


•’r 


The  inverse  Fourier  transform  of  a  function  g  G  Lj(Z  ^  x  R  is  given  by 


-  n^  n« 
g:  T  ^  X  R  ^  -d  C, 


g(d,s)  =  (2t)  ‘  ^  ^  s  L  g(k,x)e*'''=dx  e*'"'*. 

"1  \  2 

kGZ  ^  ^ 


jiift 

For  functions  G  G  S  ^  we  have  G  =  G. 

n^  n«  n.  n„ 

For  (k,<r)  G  Z  x  R  ^  and  (^,s)  G  T  x  R  we  define  the  dot  product  in 
the  natural  way,  i.e. 

"1  "2 

(k,<r).(j{,s)  =  E  k-f  +  S  s-<r-. 

i=l  ^  ^  j=l  ^  ^ 

As  to  the  set  of  points  (<J,s),  at  which  functions  are  sampled,  we  restrict 
our  analysis  to  the  so-called  sampling  lattices: 

Definition  2.1 

We  call  a  non- singular  (nj+n2,nj+n2)- matrix  V  feasible  [9]  for  sampling  on 

"1  "2  i  n  i 

T  X  R  ,  if  2Te  G  WZ  for  i  =  l,...,np  where  the  e  are  canonical  unit 

^1^^2  i 

vectors  of  R  ,  i.e.  e^  =  ^. ..  The  set 

ii-j  Rrt 

L  =  (WZ  ^  n  (T  ^  X  R 

n.+Hn 

=  {Wk|k  G  Z  ^  ^(Wk).  G  T,  i  =  l,...,nj} 
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is  called  the  sampling  lattice  generated  by  V.  The  set 

L>  =  2i\rh  ^  ^  =  {2xV^klk  6  Z  ^ 
is  called  the  dual  lattice  with  respect  to  L. 

The  simplest  example  for  a  sampling  lattice  is  the  so-called  standard 
lattice  on  T  X  R.  It  is  generated  by  a  diagonal  (2,2)- matrix  with  diagonal 
(2T/P,h)  where  P  G  N,  h  >  0  and  is  equal  to 


Ls:=  {  ^  G  T  X  R:  ^  s  =  h£;  G  Z,  0  <  j  <  pj. 


(9) 


The  condition  for  feasibility  means  that  WZ  is  2t- periodic  in  the 

first  n^  variables  and  implies  that  the  dual  lattice  is  a  subset  of  Z  ^  x  R  ^ 
which  is  the  domain  of  the  Fourier  transform  of  functions  defined  on 


^2 

T  ^  X  R 


A  given  sampling  lattice  does  not  uniquely  determine  the  generating 
W.  It  does  determine,  however,  |det  W|  as  well  as  the  dual  lattice, 
dual  lattice  L'  can  be  characterized  in  terms  of  L  by 

L'  =  {u  G  z"^  X  r"2|W  e  L:u.v/(2x)  G  Z}. 

For  a  given  lattice  L  we  define  the  lattice  constant  Cj^  by 

-(n^+nJ/2 

Cj^  =  (2x)  ^  ^  |det  W| 


matrix 

The 


where  V  is  any  feasible  matrix  generating  L. 

A  lattice  must  satisfy  a  certain  requirement  for  its  density  of  sampling 
points  in  order  to  be  a  suitable  lattice  for  sampling  a  given  function  g. 

The  following  theorem  shows  that  this  density  requirement  is  determined  by  the 

size  and  shape  of  a  set  K  c  Z  ^  x  R  ^  in  which  the  Fouier  transform  g  is 
concentrated  in  the  L^- sense:  The  translates  of  K  with  respect  to  the  dual 
lattice  must  be  mutually  disjoint.  Note  that  this  requirement  of  sparsity 
for  the  points  of  the  dual  lattice  corresponds  to  a  density  requirement  for 
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the  points  of  the  lattice  itself. 
Theorea  2-2 


u-j  Hrt  iirt 

Let  L  be  a  sampling  lattice  for  T  ^  x  R  ^  and  K  c  Z  ^  x  R  ^  such  that  the 

^2 

sets  K  +  u,  u  G  L'  are  mutually  disjoint.  Let  g  G  S  ’  and 

S^gCx)  =  Cj^  S  ;rK(x-v)g(v) 
vgL 

where  is  the  characteristic  function  of  K.  Then 


sup  |(S^g-g)(<i,s)|<2(2T) 
n.  n^ 

({},s)gT  ^xR  ^ 


(nj+n2)/2 


1 


kGZ 


The  theorem  is  a  modified  version  of  the  sampling  theorem  of  Petersen  and 
Middleton  [14] .  It  shows  that  the  function  g  is  essentially  determined  by 
its  values  on  a  sampling  lattice  L,  if  |g(k,(r)|  is  small  outside  a  set  K 
satisfying  the  conditions  of  the  theorem.  The  key  for  the  proof  of  the 
theorem  as  well  as  for  proving  some  of  the  error  estimates  for  the  filtered 
backprojection  algorithm  to  be  presented  later  is  the  following  lemma: 

Lemma  2.3 

ii'i  Hrt  jiirt  n.<  iin 

Let  L  be  a  sampling  lattice  for  T  x  R  G  g  S  ^  F  G  Lj(Z  ^  x  R 

and  F  =  (F)"'  Then  for  every  y  g  T  ^  x  R  ^ 

Cj^  S  F(y-v)G(v)  = 


*  (n.<+nn)/2  r  ■  /i  \ 

=  (2t)  1  2  j;  LF(k,(r)  S  G((k,(7)-u)e"(^’‘')-yd<T.  (11) 

n^  J/2  ugL' 
kGZ  1  ^ 

Proof:  Replace  F(y-v)  on  the  lefthand  side  of  (11)  by 
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n. 


'r  2 

keZ  ^ 

and  use  the  Poisson  summation  formula: 


S  ^  6((k,<r)-u) 


v€L 


uGL' 


(12) 


which  is  proved  in  [9] . 

Proof  of  Theorem  2.2:  Using  the  lemma  with  F  =  we  obtain’ 

l(Si,g-g)(^,s)| 


(2t) 


-  (ni-^n2)/2 


kGZ 


4  R 


S  g((k,,7)-u)-i(k,„)le'''-«'e'''-"d,. 
ugL'  J 


S  g((k,(r)-u)e^^'^e^‘''®d(r 
(k,<r)  G  ”K  uGL'  jUifO 


-  (2t) 


-(nj+n2)/2 


n^  no 


(k,(r)  G  Z  ^  X  R  V. 


-(n.+n„)/2 

<  2(2')  ^  ^  S 


lg(k,(r)|d(r 


“1  “9 

(k,(7)  G  Z  ^  X  R  ^\K 


where  we  have  used  the  fact  that  the  sets  K  +  u,u  G  L'  are  '.utually  disjoint, 
and  therefore 

S  I  S  |g((k,(r)-u)|d<r  <  S  |g(k,(r)  |d(r. 

(k,(r)  G  K  uGL',U7tO  z"^  X  r"2  \  K 

□ 
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In  the  following  we  apply  Theorem  2.2  to  the  sampling  of  the  Radon 
transform  on  T  x  R.  In  practice  one  wants  to  measure  for  a 

collection  of  different  'views',  where  each  view  consists  of  a  collection  of 
measurements  with  ^  fixed.  Since  this  corresponds  to  measuring  the  line 
integrals  for  sets  of  parallel  lines,  this  sampling  geometry  is  called 
parallel- beam  geometry.  Suitable  sampling  lattices  for  this  geometry  should 
contain  more  than  one  element  {^,s)  for  each  occurring  value  of  We  call 
lattices  with  this  property,  i.e.  with  the  property 

((iS,s)  G  L  =>  Es'  f  s:  (fi5,s')  G  L 

admissible  sampling  lattices.  The  following  lemma  characterizes  these 
lattices.  It  turns  out  that  only  a  finite  number  LM  of  different  values  of 
occurs  and  that  for  fixed  ^  we  get  a  set  of  equidistant  values  of  s,  where  the 
distance  d  between  two  nearest  neighbours  does  not  depend  on  (j). 

Lemma  2.4. 

Let  L  be  an  admissible  sampling  lattice  for  T  x  R.  Then  there  exists 
d  >  0  and  L,M,N  G  N  with  0  <  N  <  M  and  gcd(M,N)  =  1  such  that  the  matrix 


W(d,L,M,N)  = 

generates  L.  Hence  L  is  equal  to  the  set 
L(d,L,M,N) 


’2t  27rNl 
T  "LM 

.  0  d/M. 


G  T  X  R 


2irj 

^  "LM’®j^  ^  d(^+<5j/M) 


with  ^  N 


such  that  (N^j-j)/M  G  Z;  j  =  0,...,LM-1;  £  G  Z|.  (13) 

Proof:  Consider  an  integer  2x2  matrix  U  with  |det  U|  =  1.  Then  UZ^  =  Z^ 

and  the  lattices  generated  by  a  matrix  W  and  by  WU  coincide.  If  V  is 

T 

feasible  there  exists  an  integer  vector  (LpL2)  such  that 
T  T 

W(LpL2)  =  (2jr,0)  .  There  exists  an  integer  matrix  U  with  |det  U|  =  1  such 
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i 


\  T  T 

I  that  U(L,0)  =  (LpL2)  where  L  =  gcd(LpL2).  Hence  the  first  column  of  WU 

T 

•  is  equal  to  (2t/L,0)  .  Since  L  is  admissible  the  ratio  of  the  elements  of 

I  the  first  row  of  WU  is  a  rational  number.  Hence  (WL%«  can  be  written  as 

I  jj  with  gcd(M,N)  =  1.  Changing  N  by  a  multiple  of  M  does  not  chance  VUZ^. 

'  Hence  N  can  be  chosen  such  that  0  <  N  <  iMj-l.  Since  replacing  N  by  -N  gives 
the  same  new  lattice  as  replacement  of  N  by  M  -  N,  it  is  readily  seen  that 
f  allowing  for  negative  values  of  d,L,M  does  not  yield  further  lattices.  □ 


In  order  to  apply  Theorem  2.2  to  the  Radon  transform  Rf  of  a  function  f  we 
have  to  determine  a  suitable  set  K  such  that  the  righthandside  of  (10)  is 
small.  This  has  been  done  by  Lindgren  and  Ratbey  [11,  16]  who  obtained  the 
set 

=  {(k,<^)  G  Z  X  R  I  |<7-|  <  b,|kl  <  -  max(|<rl,(l-r)b)}  (14) 

which  is  shown  in  Figure  1. 
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Here  b  plays  the  role  of  an  essential  band-  width  of  f  in  the  sense  that 

is  sufficiently  small.  r  is  an  auxiliary  parameter  between  0  and  1  and 
usually  chosen  very  close  to  1.  The  following  rigorous  estimate  was  given  by 
Natterer  [12,  p.  73] : 

Theorem  2.5 

Let  f  G  C^(fi) .  Then 

E  |M(k,<r)|d<r  _<  —  ep(f,b)+||f||L^J7(r,b) 

(k,<T)  G  Z  X  R  \  KJr,b) 
where  J7(r,b)  satisfies  an  estimate 

0  <  j?(r,b)  <  C(r)e‘ 
with  constants  C(r),A(r)  >  0. 

According  to  the  last  two  theorems  the  procedure  of  determining 
appropriate  sampling  lattices  is  as  follows:  Use  a-priori  information  to 
determine  a  value  of  b  for  which  ^^(fjb)  is  small.  Choose  t  close  to  one 
(e.g.  0.99)  and  look  for  lattices  L  such  that  the  sets  K^(7-,b)  +  u,  u  G  L'  are 
mutually  disjoint. 

If  additional  information  about  f  is  available  it  may  be  possible  to 
replace  KQ(r,b)  by  a  smaller  set.  For  example  if  f  is  radially  symmetric 
Rf(k,(r)  =  0  for  k  it  0,  and  KQ(r,b)  can  be  replaced  by  the  much  smaller  set 
K(b)  =  {(0,(r):|(r|  <  b}.  This  of  course  greatly  reduces  the  necessary  amount 
of  measurements.  In  the  following  we  will  assume  that  the  only  available 
a-priori  information  about  f  is  the  essential  band- width  b  and  therefore 
always  work  with  KQ(r,b). 

For  example,  the  standard  lattice  L  given  in  (9)  is  identical  with 

O 

L(h,P,l,0)  of  Lemma  2.4.  The  dual  lattice  L'  is  therefore  generated  by  the 

b 
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mtrix 


Tfcc  sets  K  (Tj1»)  «  k,  ■  e  L”  arc  mtaally  disjoiat  iff  T  >  2i»/r  aod  li  <  x/E». 
Siace  tfiic  ladoa  traasfora  satisfies  tiic  relatioa 

=  lf(»,s)  (15) 

it  is  advaatagcotBS  to  clioosc  P  cvea,  say  P  =  2p,  p  £  R.  Jkct 
*f(2»(j^P)/(2p)2fcf)  =  If(2xj/(2p),-liO,  so  tk&t  VC  oily  liavc  to  saa^lc  for 
®  <  j  <  P- 

Applyin'  Hicorcas  2.2  and  2.>5  yields  tbc  vcll-knovn  samplin'  conditions 
p  >  b/r  and  h  <  x/b  for  the  standard  lattice: 

Theorca  2.6 

Let  f  £  C?(i!)  and  L  (h.p)  :=  L(h,2p,l,0)  be  the  standard  lattice  vitb 
p  >  b/r  and  h  <  x/b.  Then 

The  standard  lattice  is  not  the  aost  efficient  lattice  suitable  for 

saapling  Rf,  Lindgren  and  Rattey  [11,16]  found  that  a  saapling  lattice  first 

suggested  by  Coraack  [2]  which  needs  only  half  as  aany  saaples  satisfies  the 

conditions  of  the  theoreas.  This  so-called  'interlaced  lattice'  contains  the 

f2xj  -j 

saapling  points  =  1^-^,  2h^+h5jJ,  j  =  0,...,2p-l,  C  e  7,  where 


<5j  =  j 


mod  2 


"  {o 


j  odd 
j  even 


and  is  equal  to  L(2h,p,2,l)  =:  Lj(h,p).  We  obtain  the  following  conditions 
on  p  and  h: 

Theorca  2.7 

Let  f  6  C^(fl)  and  Lj(h,p)  be  the  sampling  lattice  L(2h,p,2,l)  defined  in 
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1 


(13). 


or 


taca 


If  citlier 

r/(2l»)  <  li  <  tfh  saA  p  >  (!*«Bax(x/ii-b,(l-T)b))/r 

li  <  */(2b)  aad  p  >  b/r. 

S 


(16) 

(17) 


,  -  ^^(f,b).i%  ^r,b). 

(f.S)£IxK  1  1 

« 

Proof:  Tbe  sets  K^(r.b)  -f-  a.  o  6  Lj  arc  ■utnally  disjoint  if  acd  oaly  if  the 
above  conditions  on  p  and  b  are  satisfied  Ibe  esticatc  follovs  again  frea 
Tbeorefis  2.2  and  2. -5. 

□ 

Ve  see  that  the  interlaced  lattice  has  a  spacing  of  2h  betveen  adjacent 
sasples  in  the  s- variable  and  needs  only  a  slightly  increased  value  of  p.  If 
p  is  even  the  values  Rf(dj,Sj^)  for  j  >  p  arc  redundant  because  of  (15).  In 
this  case  the  interlaced  lattice  is  nearly  twice  as  efficient  as  the  standard 
lattice.  If  p  is  O'Ld  the  data  for  j  >  p  are  not  redundant  and  we  obtain  the 
sane  data  as  from  the  standard  lattice. 

While  the  theorem  shows  that  Rf  is  determined  by  its  values  on  the 
interlaced  lattice,  it  is  not  immediately  clear,  how  good  reconstructions  of  f 
can  be  obtained  from  these  data.  The  remainder  of  this  paper  is  devoted  to 
this  question. 

3  TUB  FILTERED  BAGKPROJECTION  ALGORITHM 

In  this  section  we  describe  one  of  the  standard  reconstruction  methods,  the 
filtered  back-projection  algorithm  [15]  and  its  discrete  implementation  for 
data  sampled  on  an  admissible  sampling  lattice.  The  algorithm  can  be 
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rc^rdcd  as  a  ccmpnter  implcmcnt^tiom  of  cq,  (S)  vhicli  read 

lfjj*f(x)  =  Vjj(x.^s)Kf(#,s)dsd?,  9  =  (cos  ^,sin  d)^. 

So  the  goal  is  to  obtain  an  accurate  reconstmetion  of  instead  of  f 
itself.  The  ent-off  frequency  b  controls  the  accuracy  with  which 
approxinates  f .  If  it  is  greater  or  cqnal  to  the  essential  band-width  of 

f-lfj^*f  will  be  veiy  close  to  f .  Therefore  we  usually  choose  the  cut-off 
frequency’  cqnal  to  the  essential  band- width.  There  is  also  a  direct 

correspondence  between  the  resolution  of  the  iaage  represented  by  Vr^^f  and  the 
cut-off  frequency  b:  the  saallcst  details  discernible  in  the  inage  have  a 
size  of  roughly  2r/b- 

Sanpling  Rf  on  an  adnissible  sanpling  lattice  L  =  L(d,L,M,X)  pereits  to 

discretize  the  integrals  by  ncans  of  the  trapezoidal  rule: 

2xd  LM-1 

V*W  =:  -U  (IS) 

T  2xj 

Sj^  =  d(€+dj/M),  0^  =  (cos  ^j,sin  =  -jjji 

where  we  have  used  (13).  A  computer  implementation  based  on  (IS)  demands  the 
computation  of  the  discretized  convolution  integral 


for  all  points  x  where  the  image  is  computed  and  all  j,0  <  j  <  LM  -  1.  Since 

-this  would  be  computationally  too  demanding,  we  compute  only  Qj((Jj,llk),  with 
1  1 

II  >  0,  k  G  Z,  -  jj  <  k  <  j|  and  obtain  an  approximation  for 
linear  interpolation  with  stepsize  II.  This  means  that  we  approximate 
by 

IllQjClij.x.Oj)  :=  E  B|,(x.e.-IIk)qj{(!j,llk) 

Ic 

where  Bjj  is  given  by 


*^11^®^  -  I  0  otherwise  . 
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(19) 


In  the  folloving  vc  vill  always  assnae  that  H  is  chosen  snch  that 

d 

x/b  ^  H  =  ju  with  ■  e  M- 

A  short  coapntation  shows  that  because  of  d/(IDi)  £  N 

This  aeacs,  that  the  effect  of  the  interpolation  can  be  expressed  by  replacing 
W|^  with  the  interpolated  kernel 

k€Z 

Its  Fourier  transfon  is  given  by 

2x£ 

(WM  =  -  -it] 


sin  s  X 

where  sinc(s)  =  — - — -  Hence  IjjWjj  is  not  bandlinited.  For  \ff\  <  jj  we  have 

(Ill«b)‘(<^)  =  sinc^(H(r/2)wjj((r) 

since  we  assumed  that  b  <  x/H-  Therefore  we  can  split  IjjWj^  in  a  bandlimited 
part  Wj^Q  with  band- width  b  and  Fourier  transform  Wj^Q((r)  =  sine  (H(r/2)wjj((r), 

and  a  high-frequency  part  Wjjj  with 

2  -  r  2xfi  f  0  |(r|  <  x/ll 

“lllW  "  “bi'  -  tJ  =■  i(Irt)-(.)  k|  >  ./II. 

This  means  that  the  algorithm  computes  the  function 

2xd  LM- 1 


2xd 


LM-1 


(20) 


4  EIWOK  ESTIMATES 

From  equation  (20)  we  see  that  the  filtered  backprojection  algorithm  computes 
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a  fnaction  f^Cx)  vhich  approxiaatcs  V^^(x)  and  can  be  dccoaposcd  as 

=  fio(x)  +  fjjj(x)  where 

2rd  lM-1 
^10 W  =  Ti  ^ 

and  fm(x)  is  obtained  b\’  replacing  w^^  faj*  Wjjj.  Ve  will  show  that  f^Q 
coincides  with  St.*i  up  to  a  saall  error  and  that  K||j(x) |  can  be  kept  snail  by 
a  proper  inplenentation  of  the  algorithn.  Fron  Theoren  1.2  we  see  that 

-fo  «Lo('’«-^s)Rf(p,s)ds  dp  = 
where  h^oCO  =  2v^Ur^«Lo(i^l)- 

VloW  =  (VWb)M  (22) 

with  the  additional  low- pass  filter  Cjj  given  by 

C„(0  -  I(2')‘'siDC^(H|f|/2)  |f|<b 

"  I  0  if|>b 

compare  [12,  Theorem  V.1.2].  Hence  T^qCx)  can  be  regarded  as  approximattion 
for  Wj^Q*f(x)  =  G|j*Wjj*f  (x) .  The  discrepancy  Gjj*Wjj*f (x)-Vjj*f (x)  can  be  made 
arbitrarily  small  by  choosing  a  sufficiently  small  H. 

We  can  use  the  following  slightly  modified  result  of  Kruse  [9,  Theorem 
6.1]  to  estimate  the  error  |fj^Q(.^)-Gjj*Wjj*f(x)| : 

Theorem  4.1 

Let  f  e  C®(ft) ,g((},s)  :=  Rf({},s),  Wj^jWj^  as  in  Theorem  1.2,  L  =  L(d,L,M,N) 
an  admissible  sampling  lattice  for  T  x  R  as  given  in  (13),  and  K  c  Z  x  R  such 
that  the  sets  K  +  u,u  G  L'  are  mutually  disjoint.  Then 

2rd  LM-1 

where 
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bd  r 

-4riji  ^Js(v)l  S  I  |Jlj(«^|x|)|d(r 

^  vef.  keZ  J(k,ir)j?K,|a|<b 


with  Jj.  the  Bessel  function  of  the  first  kind  of  order  k. 

Applying  Theoren  4.1  with  Vj^  and  Wj^  replaced  by  G|j*Vjj  and  Wj^, 
respectively,  yields  the  desired  estimate  for  fj^g: 

Corollary  4.2 

Let  f,  L,K  as  in  Theorem  4.1,  f^Q  as  in  (21)  and  Vj^q  as  in  (22).  Then 

where  e^  and  62  satisfy  the  estimates  of  Theorem  4.1. 

The  sampling  conditions  of  Theorem  4.1  and  the  estimate  for  Cj  correspond 
directly  to  the  hypothesis  and  estimate  of  the  sampling  Theorem  2.2.  This 
indicates  that  the  filtered  backprojection  algorithm  is  ideally  suited  to 
e.xploit  the  advantages  of  minimal  sampling  lattices.  The  proof  of  Theorem 
4.1  given  below  employs  Lemma  2.3  in  a  crucial  way,  as  was  done  in  the  proof 
of  the  sampling  theorem. 

Proof  of  Theorem  4.1:  For  x  G  ft  define  q^:  T  x  R  R, 

=  wj.[x.[j9S  +s].  We  have 

=  S  [  q  (k,(r)g(k,(r)d<r  (24) 

keZ  JR 

and  (see  [9,  Lemma  6.1]) 

<  (4!r)-'b|J|,(Hx|)| 
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(25) 

(26) 
(27) 


Obser\’ing  that  Cj^  =  d/(lM)  the  lefthand-side  of  (23)  can  be  written  as 

SiCl  S  qx(-v)g(v). 

Applying  LcBuaa  2.3  gives  together  with  (24); 

|e(x)j  :=  |2jCj^  E  q^(-v)g(v)-Vj^*f(x)| 

1  ^  veL  ”  1 

=  I  S  f  q  (k,<r)  S  g((k,<r)-u)d(r  .  (28) 

IkeZJR  ^  uGL',u#0 

Let  i'jj  denote  the  characteristic  function  of  K.  Since  the  sets  K+u, 
u  6  L  are  mutually  disjoint,  we  obtain 
|e(x)| 


sup  „(lqx(*^>‘^)l),S  [  |g(k,<r)|d(r 


(k,<r)6K 


kez  Jr  ^  ^  vgL 

where  we  have  used  the  Poisson  summation  formula  (12)  for  the  second  term. 


Hence 


Ie(x)I  <2  sup  (Iq  (k,(r)l)  E  [  |g(k,(r)ld(r 
(k,(r)GZxR  ^  l(z7.Ja^\dV 


keZ  J(k,(r)j'K 


+  S  ls(v)l  S  [  |q  (k,<7)|d(r. 
vgL  kGZJ(k,ff)^K 

Now  the  estimate  for  e^  is  obtained  by  applying  (27)  to  the  first  term  on  the 
righthandside.  The  estimate  for  e2  results  from  obsei'ving  that  q^(k,(r)  =  0 
for  I  (7 1  >  b  and  using  (26)  for  the  second  term. 


The  results  derived  so  far  provide  the  following  procedure  for  choosing  an 
appropriate  sampling  lattice  and  obtain  accurate  reconstructions:  First  use 


a-priori  information  about  f  to  determine  the  essential  bandwidth  of  f  which 
should  also  be  chosen  as  cut-off  frequencj’  for  the  filter  If  the  cut-off 
frequency  is  too  small,  will  not  be  a  good  approximation  for  f  and  if  it 
is  bigger  than  the  essential  bandwidth,  the  reconstruction  error  02  will  be 
large.  Otherwise  62  is  not  very  critical,  as  we  will  see  below. 

Then  choose  the  sampling  lattice  L  so  that  the  translated  sets  K+u,  u  6  L' 
with  K  p  KQ(r3b)  are  mutually  disjoint-  According  to  the  results  of  section 
2  this  guarantees  that  Rf  is  properly  sampled  and  leads  according  to  Corollary 
4.2  to  a  small  reconstruction  error  Cj. 

The  application  of  Corollary  4.2  to  the  standard  and  interlaced  lattices 
is  as  follows:  As  we  have  seen  in  Theorems  2.6  and  2.7  the  required 
disjointness  of  the  sets  KQ(r,b)  +  u,  u  G  L'  translates  into  the  sampling 
conditions  p  >  b/r,  li  <  r/b  for  the  standard  lattice  and  into  the  conditions 
(16),  (17)  for  the  interlaced  lattice.  Hence  the  parameters  p  and  h  have  to 
be  chosen  accordingly.  For  the  standard  lattice  the  disjointness  of  the  sets 
Ko(r,b)  +  u,  u  G  Lg  implies  that  the  sets  Kj(r,b)  +  u,  u  G  with  the 
rectangle 

Kj(r,b)  =  {(k,(r)  G  Z  x  R:  |k|  <  b/7,|(T|  <  b} 
are  also  mutually  disjoint.  Therefore  we  can  take  Kj(r,b)  for  the  set  K  in 
Corollary  4.2.  In  case  of  the  interlaced  lattice  we  have  to  choose 
K  =  KQ(7,b)  to  achieve  optimally  sparse  sampling.  This  is  of  little 
consequence  for  e^  but  increases  e2.  To  see  this  we  reformulate  the  estimate 
for  62  and  bring  out  its  dependence  on  |x|:  Let  be  as  above, 

:=  (l-T)b/r,  /Cj  :=  b/r  and  0  :=  g(l- (r|x|)^)^/^.  Using  the  estimates  of 
[12,  p.  65]  we  obtain  the  estimate  (i  =  0,1): 
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According  to  this  estimate  we  expect  for  r  close  to  1  a  strong  growth  of 
the  error  near  the  boundary  of  fi.  Since  the  error  should  be  more 

pronounced  in  the  case  of  interlaced  sampling.  A  comparison  of  the 
reconstructions  shown  in  Figures  3  and  6  of  the  next  section  shows  that  these 
expectations  are  justified.  With  increasing  b  the  error  decays 
exponentially. 

The  performance  of  the  filtered  backprojection  algorithm  using  the 
standard  lattice  has  already  been  studied  in  [10,  12,  13,  20]  There  the 
errors  due  to  the  discrete  convolution  and  the  numerical  integration  with 
respect  to  ^  have  been  treated  separately.  It  turns  out  that  for  the 
standard  lattice  both  numerical  integrations  yield  accurate  results  when  the 
sampling  conditions  are  satisfied.  This  is  no  longer  the  case  when  the 
interlaced  lattice  is  used.  Then  the  discrete  convolutions  are  highly 
inaccurate  because  of  the  stepsize  2h  being  too  large.  The  results  derived 
in  [5]  show  that  these  errors  cancel  out  during  the  subsequent  numerical 
integration  over  T.  But  these  cancellations  may  be  disturbed  by  the 
interpolation  step  inserted  in  between  the  two  integrations.  This  provides 
an  intuitive  explanation  for  the  fact  that,  as  we  will  see  below,  the 
interpolation  is  harmless  for  standard  sampling  but  critical  for  interlaced 
sampling. 

Having  seen  that  fj^Q  represents  the  picture  we  want  to  have  it  remains  to 
clarify  under  which  conditions  |fjjj(x)|  will  be  small.  We  obtain 
Theorem  4.3 

Let  f  e  as  in  Theorem  1.2,  and  L(d,L,M,N)  be  an  admissible 


sampling  lattice.  Then 

1  rb  «  -ff  2xC]  ] 

|f„x(x)  |<  «  sup  sin'^(Ho-/2)|<r;}((r/b)|  S  f  a — j- 0  d<r. 

^  ^£T  J-b  £gZ  LL  J  J 

Proof:  With  Sj^  =  d(£+5j/M)  (see  (13))  we  obtain 


“  LM  j=0  ^eZ 


d  LM-1  r  *  i<rx.^.  iffd(f+5-/M) 

-  S  WjitW®  ^  S  e  ^  Rf (^}•,d(£+5•/M))d(r 

LM  i=n  Jr  ta  J  J 


2t  LM-1 


i<rx.0- 


LJ  JjiH- 1  r  -  -nr  2t^-^  -  i2T£<5./M 

=  IS  jfo  Jr  h’"  ■  t]'= 

where  we  have  used  Poisson's  summation  formula  (12)  with  n^  =  0,n2  =  1  and 
L  =  {d^,£  G  Z}.  With  Theorem  1.1  and  the  assumptions  (19)  we  obtain  the 
estimate 
l%I  W I 

<  (2t)^/^  sup  |wiit(<^)I  S  f  0^  -  "X"  ^ 

^SGT  Jr  M  ”  J  J 

=  sup  f  sinc^{!k/2)  «,  [j  -  ^1  S  f  ffu  -  ^1«1  in 

«Jr  ho  H  "JfcZ  It  '•i  J 


q/0  fb  sin^(iI(r/2) 

=  (2r)^/2  sup  S 

<5gT  k?tO  J-b  +  Tk  2 
2 

1  rb  „ 


)  -  ^rr  2T(kMm-^)1  T 

2  "bW.?  I  — a — ^ 


i  pu  n  -  4TC 

<  n  sup  sin^(H(r/2)  |(r?5((r/b)  I  S  f  (7  -  -j-  6  da. 
"'^iGTJ-b  Pp7.  « 


The  interpretation  of  this  estimate  is  more  involved  than  the  one  of 
Corollary  4.2.  As  an  example  we  discuss  the  case  of  standard  sampling  with 


d  =  h,  and  f  essentially  bandliinited  with  bandwidth  b  =  ir/h,  which  is  also  the 
cut-off  frequency  for  the  filter  Wj^.  Hence  2x/d  =  2b  and  therefore  the  terms 


1  «  C  r  2x^-1 

K[^  -  T-] 

"I 

with  £  ^  0  in  "  ~^J^J  |  negligibly  small.  The  term  with  i  =  0 

is  equal  to  |f(<r^)l.  In  most  applications  we  have  f(x)  >  0.  This  means 
that  |f(OI  has  a  sharply  peaked  maximum  at  the  origin.  For  |^|  close  to  the 
cut-off  frequency  b,  |f(OI  is  usually  very  small.  In  such  a  case  the 
integral  in  (30)  will  be  small  since  sin^(H(T/2)  is  small  for  \a\  «  1/H.  So 

with  standard  sampling  it  is  usually  safe  to  choose  11  =  h.  An  example  for 
this  is  the  Shepp- Logan  phantom  which  we  use  for  our  numerical  simulations  in 
the  next  section. 

If  on  the  other  hand  if(OI  is  not  small  for  |(rl  close  to  b,  the  interpolation 
stepsize  11  has  to  be  chosen  considerably  smaller  than  x/b.  Then  sin  (llu'/2) 
is  small  for  all  values  -b  <  tr  <  b.  As  an  example  for  this  case  we  will 
perform  numerical  tests  with  f(x)  =  Jj(b|x-XQ|)/|x-XQ| .  The  Fourier 
transform  of  this  function  is  constant  for  |C|  <  b  and  vanishes  for  |^|  >  b. 


5  NUMERICAL  RESULTS 

In  this  section  we  present  numerical  tests  for  the  theory  derived  so  far.  We 
will  concentrate  on  the  standard  and  the  interlaced  lattices  and  will  see  that 
the  numerical  results  can  be  understood  in  detail  using  the  theorems  of  the 
last  section.  The  theoretical  results  will  in  particular  enable  us  to  remove 
the  numerical  difficulties  with  reconstructions  from  the  interlaced  lattice 
reported  in  [9] . 

Figure  2  shows  the  first  object  we  used  for  our  tests.  It  is  a 
mathematical  phantom  due  to  Shepp  and  Logan  [18]  and  simulal^es  a  cross-section 
of  a  human  head.  Here  f  is  given  by  a  linear  combination  of  characteristic 
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functions  of  ellipses.  The  density  values  of  the  ellipses  are  the  values  of 
Rowland  [17]  multiplied  by  900  and  were  also  used  by  Kruse  [9] .  The 
displayed  values  are  those  between  1  (black)  and  75  (white).  The  biggest 
value  occuring  in  the  picture  is  900.  All  pictures  are  displayed  on  a  256  by 
256  grid.  As  essential  bandwidth  we  choose  b  =  128t  2;  402. 

Figure  3  shows  a  reconstruction  made  using  the  standard  sampling  lattice 
L  (h,p)  with  parameters 

h  =  r/b  =  1/128  p  =  404  H  =  h 
;}(<r)  =  sinc(<rT/2)jtj-.j^jj(<r)  (31) 

where  denotes  the  characteristic  finction  of  the  interval  [-1,1]- 

The  sampling  conditions  of  Theorem  2.6  are  satisfied  with  r  =  b/p  ^  0.995. 

From  the  discussions  following  Corollary  4.2  and  Theorem  4.3  we  expect  a  good 
reconstruction.  This  is  obviously  the  case. 

In  Figure  4  we  set  II  =  h/2  which  introduces  a  high-frequency  pattern  in 
the  interior  of  the  object.  According  to  Theorem  4.3  a  smaller  value  of  11 
should  lead  to  a  smaller  error  fjjj,  hence  to  a  more  accurate  reconstruction  of 
Wjj*f.  This  is  indeed  the  case.  The  high-frequency  pattern  stems  from  the 
jump  discontinuity  of  the  filter  function  f  which  causes  a  jump  discontinuity 
of  Vjj(0  1^1  =  b.  Therefore  for  this  choice  of  the  high-frequency 
pattern,  though  undesired,  is  a  true  feature  of  *  f.  It  can  be  shown  that 
if  11  =  h  =  r/b  as  in  the  previous  picture,  the  error  term  fjjj  removes  the 
discontinuity  in  the  Fourier  transform  of  the  reconstructed  function  and  thus 
causes  the  high-frequency  pattern  to  disappear. 

If  we  reconstruct  using  the  interlaced  lattice  Lj(h,p)  =  L(2h,p,2,l)  with 
the  same  parameters  (31)  as  in  Figure  3,  the  sampling  condition  (16)  is 
satisfied  with  r  =  2b/(p+b)  2:  0.998.  But  now  the  error  fjjj  caused  by  the 
interpolation  becomes  critical.  The  estimate  of  Theorem  4.3  reads 
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SiDCC  lf(^)|  is  smll  for  f^f  >  b  aod  l«c£ieso  of  &  <  %fh^  omlw  t&c  toras  wilfc 
m  <1  ia  tltc  sot  01  tkc  ri^thsiadside  are  sot  ncslisibl5f  small.  As  kc  fcave 
already  disoisscd  tfce  ter*  witli  f  =  0  will  sot  he  critical  for  as  o'i»jcct  li&e 
this  one.  For  \t[  =  1  hovcicr,  (fCC^xf/h)^!  =  |f((^£b)tf)|.  This  ateass 
that  the  integrand  os  the  ri^thaadsidc  of  (32)  might  asseme  large  valccs  for 

e% 

\a\  close  to  b,  since  sin  (ilr/2)  is  not  snail  for  these  raises  of  e  if 
H  ~  x/b.  So  we  have  to  expect  a  considerable  rcconstmctioa  error  ia  this 
case.  The  picture  of  Fignrc  -5  shows  this  error  which  nahes  the 
reconstruction  totally  useless.  The  picture  is  essentially  the  sane  as  tfce 
one  of  [9.  Fig.  8.2.  (c)].  Bat  our  estimate  (32)  also  suggests  the  folloviug 
three  ways  to  renedy  the  problca: 

1-  Choose  H  «  x/b  so  that  sin^(H<r/2)  «  1  for  \o\  <  |b|.  Vhile  this  method 
will  always  work,  the  next  two  possibilities  are  only  suitable  if  3f(OI 
is  peaked  around  the  origin. 

2.  Choose  h  snaller  than  x/b  so  that  |<7i:x/h]  stays  away  froa  0  if  \ff\  <  b. 

If  |f  I  is  peaked  around  the  origin,  even  a  smII  decrease  in  h  should  give 
considerable  inprovenent,  see  Figure  7  below.  On  the  other  hand 
decreasing  h  night  result  in  a  violation  of  the  condition  that  the 
translated  sets  K^(~,b)  u,  u  6  L'  are  nutually  disjoint  and  thus 
introduce  new  artifacts. 

3.  Choose  the  filter  such  that  lj5((r/b)|  «  1  for  \(r\  close  to  b.  Then  the 
integrand  in  (32)  remains  small  also  for  these  critical  values  of  J^rj. 
These  three  options  are  tested  in  the  following  pictures. 

In  Figure  6  we  put  H  =  x/{16b)  so  that  sin^(H(r/2)  <  0.01  for  j^I  <  b.  In 
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tB)£  imterior  t&c  giictisric  Is  vlitiEallj'  idcailicail  to  t&e  rccoastractioia  witb  tl:e 
st£!Q^2sid  lattice  aad  snail  B  s&ovai  la  Fl^rc  4.  la  partlcalar  t&c 
lilgJi-frcwscBcy  pattcois  dec  to  llic  discoatlaBlty  of  #  coincide.  Aside  frora 
t&e  wSiite  riag-type  artifact  at  tisc  boeadary  of  tbe  pictsrc  vc  laavc  a  good 
recoastnsetioa  of  lfjj*f-  Ibis  corroijorates  tbe  assertion  tlsat  tbe  large  error 
la  tbe  previous  recoastrectioj  is  caesed  by  tfec  interpolatioa,  Tbe  artifact 
at  tbe  boEffldary  can  be  traced  to  tbe  error  02  of  Corollary  4,2.  Tbe  estiatatc 
(29)  sbovs  tlsat  this  error  tern  grows  strongly  for  |x|  approaching  1.  This 
explains  that  the  artifact  appears  at  the  boundary  of  the  nnit  circle. 
Fnrthcrroorc  according  to  (29)  the  artifact  shonld  vanish  if  p  is  slightly 
iccrcascd,  corresponding  to  a  snaller  value  of  r.  Figure  S  shows  that  this 
is  indeed  the  case. 

The  second  possibility  of  iiEprovcseat.  nacely  to  decrease  h.  was  tried  in 
Figure  7.  Vc  used  again  the  paraacters  (3i)  but  set  h  =  1/150.  Hence 
Iffiz/bl  >  r/h  -  b  =  22x  if  kl  <  b  =  128x.  Hence  the  peak  of  |f (^) |  near  the 
origin  is  avoided,  and  as  expected  the  reconstruction  error  is  strongly 
reduced.  In  the  interior  of  the  picture  there  are  still  soae  disturbances. 

A  coaparison  with  the  next  picture  suggests  that  these  artifacts  result  fron 
the  overlapping  of  the  sets  Kp(r,b)  ^  u,  u  G  L'  caused  by  decreasing  h  without 
increasing  p.  In  Figure  8  we  used  again  h  =  1/150  but  increased  p  to  470,  so 
that  this  overlap  is  greatly  reduced.  In  addition  we  can  assume  a  smaller 
value  of  r  in  the  error  estimate  (29).  As  expected  from  our  reasoning  above 
the  disturbances  in  the  interior  are  removed  as  well  as  is  the  white  ring  at 
the  boundary.  The  high-frequency  structures  in  the  interior  of  the  phantom 
are  very  similar  in  Figures  7  and  8  but  clearly  different  from  the  patterns 
visible  in  Figures  4  and  6.  The  latter  ones  are,  as  we  have  seen  no 
reconstruction  errors  but  belong  to  The  difference  stems  from  the 
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intcrpolatiOB  asd  can  be  rcaored  by  rednein'  H,  c.g.  setting  H  = 

Ibis  was  essentially  done  by  Krnsc  in  [9,  Fig.  S.2.(d)3.  In  his  pictare* 
however,  the  high-frequency  pattern  appears  as  a  low-frequenQ-  disturbance 
which  was  not  conplctcly  removed  even  by  strong  smoothing,  see  [9,  Fig. 
S.2.(e)].  The  reason  for  this  is  that  the  picture  is  displayed  on  a  16C  by 
160  instead  of  a  256  by  256  grid,  which  causes  aliasing. 

The  third  method  mentioned  above,  i.c.  choosing  the  filter  v  such  that 
!FWb)|  «  1  for  \a\  close  to  b  yielded  the  best  result  for  this  phantoa. 

In  Figure  9  we  used  the  filter  ^{a)  =  cos(ffr/2),Vf_|  other 

parameters  as  in  (31) .  Ve  obtain  a  good  reconstruction  without  artifacts  and 
of  comparable  quality  as  the  reconstruction  of  Figure  3  where  the  standard 
lattice  was  used.  The  small  values  of  i5(<r/b)  for  \a\  close  to  b  seen  to 
remove  the  ring- like  artifact  at  the  boundary.  The  continuity  of  d  removes 
the  undesired  high-frequency  patterns  in  the  interior.  The  price  for  this  is 
a  slight  loss  in  resolution.  So  this  method  will  only  work  if  lf(OI  is 
already  small  for  |^1  close  to  b. 

In  summary,  the  numerical  tests  show  that  for  this  object  good 
reconstructions  using  the  interlaced  lattice  are  possible  if  the  algorithm  is 
implemented  in  the  right  way. 

Our  second  test  object  is  the  strictly  bandlimited  function 

Jl(bolx-Xol) 

=  |x->r„i . - 

T 

with  b^  =  20x  and  =  (-0.4, 0.7)  .  The  main  difference  between  the  two 
phantoms  can  be  seen  from  the  behaviour  of  their  Fourier  transforms.  The 
Fourier  transform  of  the  Shepp- Logan  phantom  is  peaked  around  the  origin  and 
relatively  small  for  |^|  close  to  b,  but  has  no  compact  support.  The  Fourier 
transform  of  the  function  above  is  constant  for  |^|  <  20t  and  vanishes  for 
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Ul  >  20r.  Since  c^(f  ,b)  =  0  for  b  >  20r  ve  should  be  able  to  obtain  very 
accurate  reconstructions  fro«  both  lattices.  Ve  will  see  that  this  is  indeed 
possible.  To  achieve  high  accuracy  we  have  to  take  into  account  that  f  has 
no  coapact  support  and  that  the  values  of 

Rf(j!,s)  =  2  sinc(20jr(s-x^.5)) 

for  |s|  >1  are  not  snail  enough  to  be  conpletcly  neglected.  For  our 
purposes  it  will  be  sufficient  to  use  the  data  for  Is)  <  6  for^  reconstructing 
f  inside  the  unit  circle.  As  neasure  for  the  accuracy  of  the  reconstructions 
we  use  a  discrete  relative  L2- error 

-jl/2r  yl/2 

ljE:=  -  %(xjj))2  E.(f(Xjj))2 

UjJ  J  U>j 

where  the  x-j  are  the  points  at  which  the  reconstruction  is  computed. 

Phantoms  of  this  kind  have  been  used  by  Rowland  [17] ,  Natterer  [12]  and 
Kruse  [9]  for  reconstructions  using  the  standard  lattice  and  by  Kruse  also  for 
reconstructing  with  the  interlaced  lattice-  All  three  authors  choose  =  0. 
This  choice  is  not  appropriate  for  our  purpose  since  it  makes  the  object 
radially  symmetric.  This  means  that  (Rf)"(k,(r)  =  0  for  k  ^  0  which  is  too 
great  a  simplification.  Rowland  demonstrated  that  it  is  necessary  to  choose 
H  «  r/b^  even  when  using  the  standard  lattice.  The  reason  for  this  can  be 
seen  from  Corollary  4.2.  The  additional  filter  G]|  is  now  critical  because 
|f(^)l  is  not  small  for  |^|  close  to  b.  While  a  choice  of  II  close  to  r/b^ 
leads  to  discrete  relative  L2-errors  around  0.4,  much  more  accurate  results 
are  possible  with  II  small.  For  our  tests  we  used  the  following  parameters: 

b  =  20t  p  =  70  II  =  T/(16b) 

M 

The  reconstruction  using  the  standard  lattice  with  h  =  r/b  =  1/20  yields  an 
excellent  result  with  L2E  =  0.0085.  Using  the  interlaced  lattice  with  the 


sa*c  para«eters  we  obtain  an  equally  good  result  with  =  0.0093.  But  in 
order  to  achieve  this  accuracy  it  is  crucial  to  satisfy  the  saapling 
condition,  i.e.  to  aake  sure  that  the  sets  K^(T,b)  +  u,  u  e  Lj,  r  =  2b/(p+b) 
are  autually  disjoint.  Increasing  h  to  1/18  increases  L2E  to  0.56.  Even 
decreasing  h  increases  the  error  as  the  following  table  shows: 


1/h 

¥ 

16 

(r.74 

18 

0.56 

20 

0.0093 

22 

0.042 

24 

0.22 

30 

0.57 

36 

0.45 

38 

0.33 

40 

0.0087 

42 

0.0072 

50 

0.0070 

The  explanation  for  this  behaviour  is  that  the  sets  Kp(r,b)  +  u,  u  G  Lj- 
are  mutually  disjoint  only  for  h  =  1/20  and  h  <  1/40.  Hence  satisfying  one 
of  the  sampling  conditions  (16),  (17)  is  absolutely  crucial  in  this  case, 
where  the  Fourier  transform  of  the  object  is  not  small  for  frequencies  close 
to  b. 
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Figure  2:  left  Figure  3:  right 


Figure  4:  left  Figure  5:  right 


Figure  6:  left  Figure  7:  right 


Figure  8;  left  Figure  9:  right 
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R.F.  MILLAR 

An  inverse  moving  boundary  problem  for 
Laplace’s  equation 

I  INTRODUCTION 

Suppose  that  one  or  more  incompressible,  viscous  fluids  flow  slowly  in  the 
narrow  gap  separating  two  plane  parallel  plates.  Such  a  configuration  is 
known  as  a  Hele-Shaw  cell.  For  steady  flow,  the  problem  is  essentially 
two-dimensional  and,  in  the  absence  of  gravitational  effects,  the  average  of 
the  fluid  velocity  across  the  stratum  is  proportional  to  the  pressure  gradient 
([11,  §330]).  Since  Darcy's  law  for  flow  through  porous  media  is  of  the  same 
form,  two-dimensional  flow  through  a  porous  medium  may  be  modelled  by  flow  in 
a  Ilele-Shaw  cell.  Because  Hele-Shaw  flow  is  more  easily  accessible  to 
observation  and  experiment  than  flow  in  a  porous  medium,  this  connection  has 
been  utilised  frequently  to  examine  phenomena  such  as  the  interfacial 
instabilities  known  to  occur  in  both;  see,  for  example,  [19],  and  the  recent 
review  articles  [1] ,  [5]  and  [18] . 

In  a  typical  Hele-Shaw  problem,  one  viscous  fluid  displaces  a  second  that 
has  different  viscosity  or  density,  and  interest  centres  on  the  evolution  of 
the  interface  between  them.  Instabilities  in  the  interface  may  arise  when 
the  displaced  fluid  is  more  viscous  than  the  other  liquid.  The  easiest 
situation  to  analyse  occurs  when  one  of  the  fluids  has  negligible  viscosity 
and  density.  In  these  circumstances,  one  may  neglect  the  motion  of  this 
fluid,  the  pressure  of  which  is  taken  to  be  zero.  We  shall  confine  attention 
to  this  most  simple  case  of  single- phase  flow  and,  although  the  relationship 
between  fluid  velocity  and  pressure  gradient  is  based  on  the  assumption  of 
steady  flow,  we  shall  assume  its  validity  in  the  time- dependent  case  as  well. 
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Ve  shall  also  employ  generally- accepted  boundary  conditions  on  the  interface, 
even  though  the  question  of  which  conditions  are  physically  correct  can  hardly 
be  considered  as  settled  ([1],  [5],  [18],  [20]). 

One  may  distinguish  between  direct  problems,  and  inverse  (synthesis  or 
control)  problems.  In  the  former,  the  real  physical  sources,  or 
singularities,  of  the  flow  are  prescribed  in  space  and  time  and,  starting  from 
some  given  initial  state,  one  wishes  to  predict  the  state  for  subsequent 
times. 

For  an  inverse  problem,  the  evolution  of  the  flow  is  prescribed,  and  the 
object  is  to  determine  what  sources  -  if,  indeed,  any  -  will  sustain  the 
motion.  Because  of  the  theme  of  this  meeting  (but  also  since  the  method  to 
be  described  is  better  suited  to  that  task  than  to  the  direct  problem),  we 
shall  consider  only  the  inverse  problem  and,  for  certain  examples,  we  shall 
discuss  whether  or  not  the  prescribed  motion  of  the  interface  can  be  realised 
in  practice. 

Conformal  mapping  techniques  have  provided  a  powerful  means  for  studying 
Hele-Shaw  problems.  For  time- dependent  (as  opposed  to  steady- state) 
problems,  the  usual  procedure  consists  in  obtaining  a  nonlinear  differential 
equation  for  the  time- evolution  of  the  unknown  analytic  function  that  maps  the 
unit  disc  onto  the  flow  region  ([1],  [6],  [7],  [8],  [9],  [15]).  It  is 
believed  that  these  methods  have  not  been  used  to  study  inverse  problems. 

In  the  present  work,  a  different  approach  is  adopted.  For  the 
single- phase  Hele-Shaw  problem,  it  will  be  seen  that  the  pressure  is  the 
solution  to  a  Cauchy  problem  for  the  Laplace  equation,  with  analytic  data 
prescribed  on  the  unknown  interface.  An  explicit  representation  exists  for 
this  solution,  analytic  in  a  neighbourhood  of  the  interface  ([12]). 

Moreover,  the  Cauchy  data  can  be  expressed  completely  in  terms  of  what  is 


called  the  Schwarz  function  S  of  the  interface  curve  ([2]).  Consequently, 
the  pressure  can  be  written  explicitly  in  terms  of  S  and  its  derivatives 
alone.  Prescription  of  S  as  a  function  of  time  is  thus  equivalent  to 
prescribing  the  evolution  of  the  interface,  and  the  singularities  required  to 
maintain  the  process  are  readily  found.  Some  of  these  will  lie  on  one  side 
of  the  interface  in  the  fluid  region;  others  may  lie  outside  the  fluid  where 
the  solution  corresponds  to  the  analytic  continuation  of  the  pressure.  (The 
Schwarz  function  has  arisen  in  previous  analyses  of  Hele-Shaw  problems;  see 
[6],  [9],  [10],  [15],  [16]  and  [17].) 

In  the  following  section,  the  problem  is  formulated.  Then  the  solution 
is  expressed  in  terms  of  the  Schwarz  function  of  the  interface  curve.  The 
effect  of  surface  tension  (T)  on  the  interface  is  included,  and  some  immediate 
consequences  are  briefly  described.  The  physical  realisability  of  solutions 
is  discussed.  Time- dependent  problems  are  examined  under  the  assumption  that 
T  =  0.  Consideration  is  given  to  interior  problems,  in  which  the  fluid 
occupies  a  bounded,  simply- connected  domain,  and  to  exterior  problems  in  which 
the  complement  of  the  fluid  domain  is  of  this  form.  Specific  attention  is 
given  to  interfaces  that  are  circles,  ellipses,  and  limacons.  The  stability 
problem  is  not  addressed.  For  the  interior  problem,  an  interface  that  is 
elliptical  for  all  time  can  be  generated  by  a  system  of  simple  sources  on  the 
interfocal  segment.  In  the  exterior  case,  an  interface  that  is  elliptical 
for  all  time  can  be  generated  only  if  the  ellipses  have  constant  eccentricity 
and  the  pressure  becomes  unbounded  logarithmically  at  infinity.  In  both 
interior  and  exterior  problems  for  limacons  there  are  singularities  in  the 
finite  plane  inside  and  outside  the  fluid  region.  For  the  interior 
problem, it  is  possible  to  generate  the  solution.  In  the  case  of  the  exterior 
problem,  however,  the  relevant  singularities  extend  to  infinity  and  it  seems 
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unlikely  that  this  solution  could  be  realised  in  practice. 

In  connection  with  these  examples,  some  consequences  of  nonvanishing 
surface  tension  are  investigated  briefly.  The  interior  solution  for  an 
ellipse  still  seems  to  be  realisable,  but  the  exterior  solution  is  not.  For 
the  limacon,  the  mathematical  solution  does  not  exist  when  T  ^  0. 

Finally,  we  mention  some  other  problems  that  can  be  formulated  in  the  same 
manner  as  those  described  here. 

2.  FORMULATION  OF  THE  PROBLEM 

Let  denote  a  domain  in  the  complex  z- plane.  Its  boundary  evolves  with 
time  t.  An  incompressible,  viscous  fluid  fills  a  region  bebween  two 
closely- spaced,  parallel  plates,  and  projects  orthogonally  onto  D^.  Suppose 
that  this  Hele-Shaw  cell  is  unbounded  and  consider  an  inverse  problem  in  which 
the  evolution  of  is  prescribed.  One  wishes  to  determine  the  sources  or 
sinks,  if  any,  that  will  produce  this  behaviour. 

It  will  be  assumed  that  is  a  simple,  closed,  analytic  curve,  oriented 
positively.  The  unit  normal  n  to  will  always  be  drawn  out  of  the  bounded 
domain  enclosed  by  C^. 

If  z  6  C^,  then  the  Schwarz  function  S  of  is  defined  by  z  =  S(z,t), 
([2]),  so  knowledge  of  S  determines  C^.  It  may  be  obtained  by  rewriting  the 
equation  f(x,y;t)  =  0  for  in  terms  of  z  and  z,  and  solving  for  z.  The 
Schwarz  function  is  analytic  near  C^. 

Let  the  angle  i)  be  defined  in  terms  of  the  unit  positive  tangent  vector  t 
to  by  t  =  (cos  i),sm  f).  Then  n  =  (sin  ^,-cos  it),  and  the  derivative 
S^(z,t)  =  z  €  Cj  [2,  (7.8)].  Ve  shall  define  (Sj,{z,t))‘''^  by 

:=  e'*>*,  zee., 

from  which  it  follows  that  dz/ds  =  S’^/^  on  C..  Thus  with  sj/^  is  associated 

2  u  Z 
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»  particalar  oricataticM  of  C^. 

i  discassioB  of  tte  lele-Skw  eqaatioos  is  £im  is  [19].  If  v  devotes 
the  flaid  velocitv  ia  thea,  if  sraritatioaal  effects  are  aegligible, 

a  =  -  ^  Ta,  (2.1) 

ia  vhich  E  :=  fijft/d  vhere  d  is  the  distaace  betveea  the  plates,  p  is  the 
coefficieat  of  viscosity  of  the  flaid,  aad  a  is  the  flaid  pressare. 
Incoopressibility  iaplies  that 


Aa  =  0,  2:  €  D^, 

(2.2) 

except  at  singularities. 

.4t  a  poiat  m  C^, 

n  =  iTjc, 

(2.3) 

and  v.n  =  V_,  so 
n* 

Ba 

(2.4) 

here  T  is  the  surface  tension  coefficient,  assnacd  to  be  son- negative  for 
iaaiscible  fluids,  k  :=  dp/ds  denotes  the  curvature  of  C^,  and  is  the 
velocit}'  of  at  z  in  the  direction  n.  The  sign  aabiguity  in  (2.3)  is 
resolved  in  the  following  Banner:  for  interior  problcas  in  vhich  contains 
the  upper  sign  is  chosen;  for  exterior  probleas  ve  choose  the  lover  sign. 
(This  sign  aabiguity  is  a  result  of  the  unique  definition  of  for  both 
interior  and  exterior  probleas.) 

3.  SOLUTION  lEPRESENTATION 

Equations  (2.2)  to  (2.4)  define  a  Cauchy  problca  for  Laplace’s  equation. 

Since  G^.  and  the  data  are  analytic,  a  unique  solution  exists  in  a 
neighbourhood  of  G^.  It  is  our  intention  to  determine  the  singularities  of 
this  solution  when  G^  is  prescribed.  The  goal  may  be  attained  by  using  a 
representation  for  the  solution. 
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Ssj^iose  tkit  D  is  a  real  aelsliiboiriiood  of  a>  a>al|1;ic  arc  C.  Let  S 
«icaote  the  Schvaiz  fuctioa  of  C,  a>d  suppose  that  ■  is  hanawic  on  D. 

Define  z  :=  x  -s-  iy,  z  :=  x  -  iy,  where  x  and  y  nay  be  coaplex.  Asscae  that 
C  is  paraaetrised  hy  ardength  s  as  z  =  z(s),  and  let  s  =  s(z)  be  the  inverse 
relationship.  If  «  =  v(s),  dikidu  =  v(s)  on  C,  then  ([12])  n(x,y)  =  U(ziz)- 
tfhcrc 


♦  1  *1  rs(z) 

U(z-z  )  =  2[v(s(z))  +  v(s(z  ))]  +  2  ^ 

^  ^s(z  ) 


here  an  overbar  denotes  the  conjugate  function  and  the  integral  is  in  the 
coaplex  plane. 


If  V(z)  :=  v(s(z))  and  h’(z)  :=  v(s(z))ds/dz.  then  (3.1)  beconcs 

U(z,z*)  =  2[V(z)  ^  V(S(z*))]  ^  5  if  ♦  V(Od(, 

■»S(z  ) 


(3.2) 


and  z  =  z  gives 


«(x,y)  =  J[V(z)  +  V(S(z)]  4  i  r  V(Od(.  (3.3) 

^  ^  JS(z) 

The  quantities  k  and  V  in  (2.3)  and  (2.4)  nay  be  e.xpressed  in  terms  of  S 


and  its  derivatives;  for,  by  [2, (7. 17)], 


/c  =  -  iS,,S, 
2  2*2/ 


3/2 


for  z  G  G^.  Also,  if  now  z(t)  :=  x(t)  +  iy(t),  where  (.x(t),y(t))  always  lies 


on  C^,  then  differentiation  of  z(t)  =  S(z(t),t)  immediately  gives  v.n,  and 


Insertion  of  the  Cauchy  data  into  (3.3)  yields 


(3.4) 
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1  1 
1  /  ZZ'  ^  ^  1  fZ 

»(x.y;t)  =  *  J  T  »e  —— — 375  -  j  ‘L _ Sj(C,t)<i<.  (3.5) 

and  the  solution  is  expressed  completely  in  terms  of  S  and  its  derivatives. 

It  may  be  shown  that  the  integral  in  (3.5)  is  real.  The  contour  of 
integration  is  obtained  by  continuously  deforming  the  path  for  a  point  near 
into  that  for  the  given  point  z,  and  will  depend  on  the  mapping  properties  of 
S.  If  S  is  multivalued,  this  integral  will  be  on  the  corresponding  Riemann 
surface.  It  is  possible,  however,  to  show  that 

f  Sj((,t)d( = r  s  ((,t)<i( + r  'sj((.t)(i(, 

in  which  z^  is  any  point  on  C^.  (A  similar  observation  has  been  made  in 

[4].)  Consequently  (3.5)  may  be  rewritten  as 

1  S„(z,t) 

{X  ZZ  rZ  ^ 

L  )i  ^0 

and  if  a  complex  potential  w  is  defined  by 

1 


X  ZZ^  ^  '  rZ 

:=  i  2  ro  f,  ,n3/2  L 


[Sj(z,t)]'  -0 

(to  which  any  purely  imaginary  function  of  t  may  be  added),  then 

u(x,y;t)  =  Re  w(z,t). 


(3.7) 


Now,  =  \  -  iUy,  so 

"x  -  i"y  =  *  (3.8) 

in  which  {Sjz}  is  the  Schwarzian  derivative  of  S.  From  (3.8)  for  z  6  C^,  and 

1 

the  result  (7.23')  of  [2]:  d/c/ds  =  2i{S,z},  together  with  (2.1)  and  3.4),  we 
deduce  that 

T  dK 

and  verify  that 
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v.n  =  (3.10) 

on  C^.  Thus,  if  T  #  0  there  is  a  tangential  component  of  fluid  velocity  on 
unless  d«/ds  =  0.  This  component  vanishes  if  is  a  circle;  but 
generally  there  is  a  discontinuity  in  tangential  velocity  across  C^,  which 
corresponds  to  a  layer  of  vorticity  on  C^. 

The  results  (3.9)  and  (3.10)  refer  to  behaviour  on  C^.  Suppose,  however, 
that  we  are  interested  in  fluid  velocity  at  some  point  P  interior  to  the 
fluid.  Then  (3.8)  is  valid,  in  which  S  refers  to  C^.  Let  C  be  any  simple, 
closed,  analytic  curve  through  P,  lying  in  the  domain  bounded  externally  by 
and  having  positive  orientation.  Suppose  that  t'  and  n'  are  the  unit  tangent 


and  normal  vectors  to  C  at  P.  Thus,  if  t'  =  (cos  sin  ^'),  then 

e^^  =  (S')'^^^,  where  S'  is  the  Schwarz  function  of  C  ;  this  result  depends 


only  on  the  direction  of  C  at  P.  Consequently,  at  P, 


v.(t'+in')  =  ±  ^(S^; 


in  which  the  result  S^^ISjz}  =  -2 — «  been  used. 

dz^ 

If  T  =  0,  the  complex  potential  and  velocity  are  determined  by  S^  alone; 
in  particular,  time  independent  singularities  of  S  do  not  play  a  role.  To  be 


more  precise,  let  z  G  D^i  for  0  <  t'<  t,  and  suppose  that  there  are  no 
singularities  of  pressure  in  D^|.  Then,  from  (3.8)  with  T  =  0, 

Iq  ■  S(z,t)), 

from  which  it  follows  that  the  singularities  in  of  S(z,t)  coincide  with 
those  of  S(z,0),  are  constant  in  time,  and  do  not  affect  w^(z,t).  These 
points  have  been  noted  previously  ([9],  [10]),  and  used  by  Ilowison  to  simplify 
the  integration  of  a  system  of  differential  equations  that  determine  the 
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evolution  of  for  an  exterior  problem.  On  the  other  hand,  if  T  0,  then 
zeros  of  S_  become  important.  For,  if  S  (z,t)  =  (z-z^)"H(z,t) ,  where  H  is 
analytic  near  z^,  H(zpt)  ^  0,  and  n  is  a  positive  integer,  then,  near  z^, 
Szz'S^^/^  ~  n(z-Zj)‘^"^/^”H(zj,t)’^^^;  thus  w  has  a  branch  point  at  z^  of 
order  at  least  -3/2  if  n  is  odd,  and  a  pole  of  order  at  least  two  if  n  is 
even.  Although  mathematically  significant,  such  singularities  are  likely  to 
be  difficult,  if  not  impossible,  to  synthesise  physically.  In  general,  then, 
such  an  S  will  not  be  admissible,  and  the  corresponding  will  not  be 
realisable.  Some  illustrative  examples  will  be  given  later. 

The  elimination  of  possible  solutions  by  surface  tension  is  consistent 
with  an  observation  in  [6].  There,  when  T  =  0,  it  was  shown  that  a  cusp 
develops  on  the  fluid  boundary  when  a  zero  of  the  derivative  of  a  function 
that  maps  the  fluid  region  conformally  and  1-1  onto  the  unit  disc  reaches  the 
unit  circle  from  its  exterior.  Such  a  zero  may  be  shown  to  correspond  to  a 
zero  of  S^,  in  the  fluid  region, 

4.  PHYSICAL  REALISABILITY  OF  SOLUTIONS 

In  studying  inverse  problems,  one  objective  is  to  determine  all  singularities 
of  the  solution.  Some  of  these  may  be  internal  singularities  that  lie  in  the 
region  occupied  by  fluid  and  others  may  be  external  singularities  that  exist 
in  the  analytic  continuation  of  the  solution  beyond  the  fluid  boundary 
G^.  For  a  given  interface  C^,  into  which  category  a  particular  singularity 
falls  depends  on  whether  we  are  considering  an  interior  problem  or  an  exterior 
problem. 

The  external  singularities  may  be  regarded  as  image  singularities  that  are 
induced  in  the  solution  and  depend  on  the  form  of  C^.  This  interpretation  is 
common  in  electrostatics,  where  the  effect  in  the  region  of  interest  of  the 


111 


image  of  a  source  in  a  conductor  is  produced  by  an  appropriate  induced  charge 
distribution  on  the  conductor.  In  the  Hele-Shaw  probem,  the  fluid  external 
to  reacts  on  against  the  internal  fluid  to  create  an  effect  equivalent 
to  that  of  the  external  singularities,  thus  resulting  in  the  required  pressure 
on,  and  motion  of,  C^.  Consequently,  we  need  not  account  explicitly  for  the 
effect  of  external  singularities,  except  to  note  that  analyticity  of  will 
be  destroyed  at  time  t  at  a  point  where  such  a  singularity  meets  ([6], 

[10],  [15]). 

A  further  aim  is  to  decide  whether  the  solution  is  uniquely  determined  by 
its  internal  singularities  and,  if  so,  whether  they  can  be  realised  in 
practice.  For  an  interior  problem,  it  was  shown  originally  by  Richardson 
([15],  [16])  that  if  the  internal  singularities  are  simple  point  sources,  then 
all  the  complex  moments  of  are  determined,  given  the  initial  configuration. 
From  these  is  found  the  part  of  S  that  is  analytic  outside  and  zero  at 
infinity.  ■  (Richardson's  analysis  would  appear  to  generalize  to  other  types 
of  singularity.)  It  was  also  pointed  out  in  [16]  that  knowledge  of  the 
moments  of  alone  does  not  necessarily  determine  uniquely,  as  shown  by 
the  example  of  Sakai  ([21]).  Nevertheless,  from  the  fact  that  functions  and 
domains  in  the  Hele-Shaw  problem  depend  continuously  on  time,  Richardson  was 
led  to  conjecture  that  some  such  domains  would  be  uniquely  determined  by  their 
moments,  and  thus  by  their  internal  singularities.  Specifically,  he  has 
determined  the  evolution  of  when  there  are  as  many  as  four  internal  simple 
sources.  His  method  involves  the  conformal  mapping  procedure;  it  becomes 
more  and  more  complicated  as  the  number  of  singularities  increases,  and 
impractical  if  the  singularities  of  S  inside  are  not  poles  ([10]).  (Note 
that  S  is  meromorphic  inside  if  and  only  if  the  mapping  function  occurring 
in  this  procedure  is  rational  ([2,  p.  158]).) 
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In  the  examples  that  follow,  we  shall  assume  that  internal  singularities 
in  the  form  of  simple  point  sources,  doublets,  or  continuous  distributions  of 
simple  sources  on  finite  curves  can  be  synthesised  physically  and,  even  when 
T  0,  like  Richardson,  we  shall  conjecture  that  the  internal  singularities 
and  initial  configuration  uniquely  determine  the  evolution  of  C^. 

5.  TIKE- DEPENDENT  INVERSE  PROBLEMS  WITH  T  =  0 
Here  we  consider  time- dependent  problems  in  which  surface  tension  is 
neglected.  Many  such  solutions  have  been  obtained:  see,  for  example,  [6], 
[S]  j  [9])  [10],  [15],  [16],  [17].  In  Richardson's  work,  solutions  are 
obtained  that  can  be  generated  by  realistic  sources.  It  seems  less  likely, 
however,  that  many  of  the  other  solutions  can  be  realized  in  practice,  even 
with  the  neglect  of  surface  tension  effects. 

Let  us  begin  with  a  simple  example.  Suppose  that  is  the  circle  of 
radius  r(t),  with  centre  at  0.  Then  S(z,t)  =  r  (t)/z,  and  (3.5)  gives 

u(x,y;t)  =  -Krr[log  z  -  log(r^/z)], 

in  which  the  principal  branch  of  the  logarithm  is  selected.  Then 
2 

log  z  -  log(r  /z)  =  21og(|z|/r),  and 

u(x,y;t)  =  -2Krr  log(|zl/r),  (5.1) 

a'result  that  follows  also  from  (3.6). 

This  solution  is  valid  whether  D^  is  bounded  internally  or  externally  by 
G^.  For  the  interior  problem,  0  <  [zj  <  r,  and  one  finds  u  ^  0  accordingly 

as  r  ^  0.  When  r  >  0,  fluid  is  injected  at  0  and  extraction  takes  place  if 

r  <  0.  The  fluid  velocity  is  radial,  the  outward  component  being  rr/|z|. 

The  source  is  at  0,  and  the  flux  across  any  curve  enclosing  0  is  equal  to  the 
time  rate  of  change  of  the  area  of  D^. 
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When  lies  outside  C^,  |z|  >  r  and  u  ^  0  accordingly  as  r  ^  0.  If 

r  >  0,  fluid  is  removed  at  infinity;  fluid  is  injected  if  r  <  0.  Again,  the 

velocity  is  radial  and  equal  to  rr/lz|,  and  there  is  a  logarithmic  source  of 
pressure  at  infinity.  In  practice,  this  flow  could  be  generated  in  the 
region  between  a  circle  centred  on  0  in  jz]  >  r  and  the  circle  \z\  =  r  by 
setting  up  the  appropriate  pressure  difference  between  these  circles,  the 
pressure  in  |z|  <  r  being  zero. 


5.1.  Interior  problem  for  ellipses 

We  consider  a  rather  more  difficult  example,  in  which  the  curves  C^(t  >  0) 

evolve  in  the  form  of  a  family  of  ellipses. 

0  0  2  2  2  2 

Let  be  the  ellipse  b  x  +  a  y  =  a  b  ,  where  a  and  b  are  prescribed 
functions  of  time,  and  a  >  b  >  0.  The  Schwarz  function  for  is 


a  +b 


2.2a/2 


S(z,t)  S  — 2"  z  -  “  (z  -c  ) 
c  c“ 


(5.2) 


2  2  2 

([2,  (5.13)]),  in  which  c  :=  a  -  b  ,  and  the  complex  plane  is  cut  from  -c  to 
c  with  (z^-c^)^/^  >  0  for  z  >  c. 

We  shall  consider  the  velocity  field  v  of  the  fluid  in  the  domain 
contained  within  G. .  It  is  given  by  (2.1)  where,  by  (3.8), 


d 


U  rcl  TU  0  0  1 

"x  ■  “y  =  3t  (— r]  (“^  -'=  ) 


d  r2abi 


1  2abc- 


(z2-c2) 


(5.3) 


Integration  of  (5.3)  on  a  curve  G  enclosing  the  segment  (-c,c)  shows  that 
the  flux  across  G  is  the  rate  of  change  of  the  area  of  D^. 


The  analytic  function  u  -  in  may  be  represented  by  its  Cauchy  integral 

X  y 

and  so  related  to  its  singularities.  Deformation  of  the  path  of  integration 
around  singularities  and  into  a  circle  at  infinity,  followed  by  some 
simplification,  leads  to  the  following  result: 


“x  ■  «y  =  at 


a-bi 


a+E 


K  fC 


r  r/7«2  A 

?  )  at 


r2abi 


2abci  d^ 

2.  A  “^J  Fi 


vic"-n 


and,  since  u„  -  iu ,  =  w^,  we  find  that 

X  y  z 


w(z,t) 


vJ 

at 


d  ra-b 


a+E 


2K 


rC  r 


-C 


)at[;:2 


abci 


log(z-Od{*g(t). 

(5.4) 


Here  any  branch  of  the  logarithm  may  be  selected,  since  only  u(=  Re  w)  is  of 
interest. 

From  (5.4)  it  is  seen  that  w  is  generated  by  internal  simple  sources  along 
the  interfocal  segment,  while  the  term  in  z  arises  from  external  sources  at 
infinity.  In  general,  c  is  dependent  on  time  and  the  internal  sources  are 
not  stationary. 

Ve  consider  two  special  cases,  but  only  one  in  detail.  In  the  first,  it 

d  i-a-  b-j 

is  assumed  that  the  ellipses  have  constant  eccentricity.  Then  =  0, 

d  rabi 


3t 


=  0,  and  (5.4)  simplifies  to 


2Kab  rC  log(z-^) 

The  function  g  may  be  chosen  so  that  u  =  0  on  C^,  and 

Wa+b  1  1  pC  log|z-{l 

-f 


(5.5) 


From  this,  one  sees  that  the  motion  is  generated  by  a  continuous  distribution 


2  2 

of  simple  sources  on  (-c,c)  with  density  -2K(ab/T)//(c  -F)>  "C  <  ^' <  c. 
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When  z  lies  on  the  segment  (-c,c),  it  is  found  that  the  integral  in  (5.5) 
1 

is  equal  to  i  log (2  c) ,  so 

.  ra+bi 

u(x,0;t)  =  2Kab  log  »  -c  <  x  <  c;  (5.6) 

here  the  logarithm  is  constant  in  time,  and  u  is  constant  for 
-c<x<c,  y  =  0,  for  a  given  time. 

Let  C  be  an  ellipse  in  that  is  confocal  with  C^.  Then  it  may  be 
shown  that  u  is  constant  on  C.  More  precisely,  if  a'  and  b'  are  its 


semi- major  and  -minor  axes,  respectively,  then 

.  <•  a+b 

u(x,y;t)  =  2Kab  logj^^r^J, 
which  reduces  to  (5.6)  if  C  collapses  onto  (-c,c) 


z  6  C, 


Since  u  is  constant  on  C,  v  is  normal  to  C.  A  calculation  based  on 
(3.11)  with  T  =  0  verifies  this  and  gives 

a'ab 

V(a'  -c  X  ) 

in  which  n'  is  the  unit  normal  to  C.  This  determines  the  analytic  function 

u  -  iu  on  C  and,  hence,  everywhere.  In  principal,  then,  if  this  velocity 

A  y 

is  impressed  across  C,  the  correct  motion  will  be  established  inside  C^;  but 
because  analytic  continuation  is  unstable,  small  errors  in  impressed  velocity 
on  C‘  will  lead  to  large  errors  elsewhere. 

For  the  second  case,  we  assume  that  the  ellipses  are  confocal.  Then,  on 
omitting  details,  it  is  found  that 

u(x,y;t)  =  2(y^-x‘^)  ^  ^  ^(ab)  /(c2-e^)log|z-e|de 


a'^b  d  rb.  1 


U  ru-j  1  r  ra+bi  l-i  d 

at  [jJ  J  -  2J  3t 


It  may  be  seen  that  u(x,0;t)  depends  on  x  when  -c  <  x  <  c,  in  contrast  to  the 


result  in  the  previous  example.  In  fact,  the  curves  on  which  u  is  constant 
are  confocal  ellipses  only  when  the  ellipses  are  members  of  a  family  with 
constant  eccentricity. 

5.2  Exterior  problem  for  ellipses 

Suppose  that  the  domain  occupied  by  fluid  is  outside  C^,  with  extraction  from, 

or  injection  into,  at  infinity.  It  is  assumed  that  S  (and,  hence,  S^)  has 

no  singularities  in  the  finite  plane  outside  C^;  thus  the  only  sources  in  the 

fluid  region  are  at  infinity.  Ve  shall  also  suppose  that  u  -  iu  behaves 

like  z'^  as  |z|  h  oo,  so  that  the  time  rate  of  change  of  the  area  inside  is 

bounded;  any  greater  growth  at  infinity  seems  to  be  unreasonable  on  physical 

grounds.  Consequently  is  of  order  z"  at  infinity.  Finally,  it  will  be 

assumed  that  S  is  at  most  of  order  z  as  |z|  oo. 

Subject  to  the  assumptions  on  S  alone,  it  may  be  shown  that  C^(t  >  0) 

determines  a  family  of  ellipses  ([22],  [23],  [14]),  for  which,  in  general, 

is  of  order  z  at  infinity  and  the  behaviour  of  u  -  iu  is  inadmissible.  (It 

X  y 

is  conjectured  that  u  -  iu  would  have  inadmissible  behaviour  at  infinity  if 

X  y 

S  grew  more  rapidly  than  z,  but  we  can  offer  no  proof.)  Only  if  the  family 

.  1 

of  ellipses  has  constant  eccentricity  will  be  of  order  z  ,  as  required, 
and  clearly  must  be  a  member  of  this  family. 

Then,  subject  to  the  assumptions  above,  one  concludes  that  can  sweep 
over  every  point  outside  if  and  only  if  C.  is  an  ellipse,  and  C.  (t  >  0)  is 
an  ellipse  with  the  same  eccentricity  as  G^.  This  result  has  been  given 
previously  in  [8] ,  and  generalised  to  K”  in  [3] . 

To  discuss  the  possibility  of  synthesising  this  solution,  the  form  of  u  at 
infinity  is  needed.  With  notation  unchanged  from  the  previous  section,  (3.7) 
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and  S^(z,t)  =  2ab(z^-c^)'^/^  give 

.  pz  d( 

w(z,t)  =  -2Kab  ..2  2x1/2’  ^ 

Jzq  (C  -c  )  / 

Since  u  is  independent  of  z^  we  may  take  z^  =  a.  Then  the  integral  is 
expanded  for  |z|  >  c  to  yield 

1  3 

w(z,t)  =  -2Kab[log  z  -  4  c^z  ^  -  ^c^z  ^  +  ...]  +  h(t), 
in  which  log  denotes  the  principal  branch,  and  h  is  real.  Thus 

u(x,y;t)  =  -2Kab  loglz|  +  h(t)  +  0(|z|  ^),  |z|  -»  00. 

If  is  a  circle,  so  c  =  0,  this  result  reduces  to  (5.1).  There  are  no 
internal  singularities  of  pressure,  but  there  is  a  logarithmic  singularity  at 
infinity.  The  external  singularities  lie  along  the  interfocal  segment  of  the 
ellipse  C^. 

In  practice,  instead  of  a  singularity  at  infinity,  a  source  or  sink  of 
pressure  would  be  impressed  along  some  closed  curve  C;  would  be  bounded 
internally  by  and  externally  by  O'.  If  C'  were  chosen  to  be  a  large 
circle  centred  on  0  then,  at  least  initially,  the  desired  flow  would  be 
generated.  Alternatively,  if  O'  were  an  ellipse,  instantaneously  confocal 
with  and  with  semi-axes  a'  and  b'  (a'  >  a,  b'  >  b),  then  the  pressure  on  G' 
is  given  by  (5.7).  The  maintenance  of  this  pressure  on  G'  would  produce  the 
flow  in  the  region  between  G'  and  but,  because  G'  changes  with  time,  this 
would  be  difficult  to  realise  physically. 

5.3.  Interior  and  exterior  problems  for  limacons 

As  a  further  example,  it  is  assumed  that  G^  is  a  limacon.  This  has  been  a 
popular  choice  in  earlier  work  that  uses  the  conformal  mapping  technique  to 
generate  solutions,  because  the  mapping  function  is  a  quadratic  polynomial; 
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but  no  attempt  seems  to  have  been  made,  to  relate  the  solutions  to  their 
singularities. 

In  polar  coordinates,  with  the  pole  of  the  system  at  z  =  -a,  has  the 

form 


r  =  1  +  2a  cos  0; 


here  a  is  a  function  of  t  and  it  is  assumed  that  0  <  a  <  25  so  that  is 
simple.  Then  ([2,  p.54]) 

(R+1) [2z+a(R+l)] 


S(z,t)  = 


4z" 


(5.9) 


in  which  R  :=  (l+4az)^/^,  and  R  >  0  for  z  >  -l/(4a).  Since  a  <  2>  the  branch 
point  lies  outside  C^,  and  the  complex  plane  is  cut  along  the  negative  real 
axis  from  -l/(4a)  to  -00.  If  a  =  0  for  t  =  0,  then  is  a  circle  of  unit 
radius,  centred  on  0.  The  branch  point  moves  to  the  right  as  a  increases, 


1 

and  meets  when  a  =  2* 

The  family  of  limacons  considered  here  is  a  one- parameter  subset  of  the 
original  two- parameter  family  considered  by  others  ([6],  [7]).  In  the 
earlier  work,  the  two  parameters  are  related  in  such  a  manner  that  the 
corresponding  has  only  a  simple  pole  at  z  =  0.  Here,  in  contrast,  will 
have  a  pole  of  order  two  at  z  =  0.  More  precisely,  for  a  >  0,  S  has  a  pole 
of  order  two  at  z  =  0  and  no  other  singularities  inside  C^.  Near  z  =  0, 

a  2a^+l  3 

S(z,t)  =  —  +  -rr-  +  (a-a  )  +  ...,  (5.10) 


and 


S^(z,t)  =  -2  +  —  +  (l-3a^)a  +  ...  .  (5.11) 

z 

By  employing  (3.8)  and  the  Cauchy  integral  representation  for  u^  -  iu^,  we 
obtain 
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(5.12) 


j-a  4aa  a  pi/ (4a)  2^^+6a^+l  ■ 

"x- ^  =  — -sj.^  R 


Thus,  to  within  an  additive  function  of  t  alone, 


and 


r  ax  .  a  pi/ (4a)  2^  +6a^+l 

u(x,y;t)  =  K  -  4aa  loglz|  -  - - 

Lx^V  J-oo  r|l+4a^| 


1/2 


logk-z|d^]. 


Suppose  now  that  the  fluid  region  is  interior  to  C^.  The  internal 
singularities  at  z  =  0  correspond  to  a  simple  source  and  a  doublet.  The 
external  singularities  are  a  distribution  of  simple  sources  on  the  negative 


1 

real  axis  from  -w  to  -l/(4a).  If  0  <  a  <  ^,  the  motion  will  be  generated  by 

1 

the  internal  sources;  when  a  =  2>  an  external  singularity  meets  C^,  which 
loses  analyticity  at  this  point  and  the  solution  breaks  down. 

If  the  fluid  is  external  to  C^,  then  the  internal  singularities  are  simple 
sources  on  (-<»,- 1/ (4a)),  and  the  external  singularities  correspond  to  a  simple 
source  and  a  doublet  at  the  origin.  Starting  from  an  admissible  initial 
configuration  (for  example,  a(0)  =  0  so  that  is  the  unit  circle),  this  flow 
could  be  generated  by  the  source  distribution  on  (-(»,-l/(4a));  in  practice, 
this  source  configuration  would  be  difficult  to  set  up. 


6.  TIME- DEPENDENT  PROBLEMS  WITH  T  ^  0 

We  shall  now  briefly  re-examine  the  examples  discussed  earlier,  but  with 
surface-tension  effects  included.  It  has  already  been  noted  that  zeros  of  S^ 
in  the  flow  region  are  not  permissible,  so  any  such  cases  will  be  omitted. 

It  goes  almost  without  saying  that  in  general  one  effect  of  surface 
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teasioa  is  to  coaplicate  tlie  analysis.  Mat  if  is  a  cJrde,  f&e  CBrsainsie 
K  is  ctwstaatx  and  tlie  coaditioa  (2.3)  aorely  raises  or  lovars  llie  pressare  m 
by  a  coBStaat  aaowt.  Tkas  sarface  teasioa  does  not  affect  tke  aotioa  ia 
this  case,  aa  observation  aade  previoasly  in  [24]. 

The  coaplex  potential  is  given  hy  (3.7),  froa  vhich  it  follows  that 

V  -  (6-1) 

Once  again  n^  -  in^,  aay  be  related  to  its  singalaritics  In*  aeans  of  the  Caachy 
integral  representation.  In  addition  to  singalarities  that  arise  vhea  T  =  0, 
the  zeros  of  nov  play  a  role. 

For  the  ellipse,  froa  (5.2)  ve  find  that  S  has  tvo  sinplc  zeros  at 
0  0  0 

z  =  ±^a"+b")/c.  Since  a(a-c)  +  b"  >  0,  it  follovs  that  these  lie  ontsidc  the 

ellipse.  Thus  the  term  in  (3.7)  or  (3.S)  that  depends  on  T  has  tvo  branch 

points  inside  and  tvo  ontsidc  C^,.  Since  S^  has  zeros  ontsidc  the  ellipse, 

for  reasons  nentioned  earlier  ve  shall  disregard  the  exterior  problca.  For 

the  interior  problea,  additional  cuts  are  nov  nade  in  the  coaplex  plane 

outside  froa  ±(a  ^-b  )/c  to  ia>.  In  addition  to  the  integral  arising  vhen 

T  =  0,  the  representation  for  w  vhen  T  #  0  vill  contain  an  integral  on  (-c,c) 

and  integrals  on  these  new  cuts.  Their  integrands  vill  behave  like 
2  2-3/2 

[z±(a  +b  )/c]  '  near  these  points,  and  the  integrals  will  involve  second 
derivatives  of  the  potentials  of  siaple  source  distributions.  Since  these 
are  external  singularities,  we  conjecture  that  the  flow  is  realisable,  in 
accordance  with  the  discussion  in  section  4. 

In  the  case  of  a  limacon,  from  (5.9)  it  is  found  that  S  has  only  one 

f  z 

2  ^ 

zero:  z  =  -2a(l-2a  );  because  0  <  a  <  2J  this  lies  inside  C^,  and  the 

corresponding  term  in  w  will  have  a  branch  point  of  order  -3/2  at  this  point. 

Since  it  is  impossible  to  define  u  -  iu  by  (6.1)  as  a  single-valued  analytic 

X  y 
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fncccticw  ia  a.  acig^lMnnieood  of  ve  CMclade  tliat  the  liaacpn  is  cot  a 
solsticw  for  cither  iatcrior  or  exterior  prohleas  if  T  #  0.  This,  thca,  is  a 
specific  exaaplc  that  illastratcs  hoc  sarfacc  teosion  liaits  the  class  of 
possible  solatioBS. 

7.  OHiCLWIllG  lEUlIS 

The  CBC- phase  Rclc-Shau  flcv  prpblcii  oo  an  nnboaeded  doaaia  caa  be  foraulatcd 
as  a  Cauchy  problea  for  the  Laplace  equation.  By  adopting  this  point  of 
viev,  and  using  an  explicit  reprcscatation  for  the  solution,  the  required 
hamonic  function  is  expressed  coaplctcly  in  tems  of  the  Schwarz  function  of 
the  interfacial  cur\'c.  This  fora  of  the  solution  is  well- suited  to  the  study 
of  inverse  probleas,  in  which  the  Schwarz  function  is  prescribed.  The 
realisability  in  practice  of  the  flow  when  the  interface  is  circular, 
elliptical,  or  in  the  for*  of  a  liaacon,  has  been  exasined  by  relating  the 
Schwarz  function  to  its  singularities,  with  or  without  the  consideration  of 
surface  tension  effects. 

Other  problcHS  can  be  studied  in  this  way.  These  include  steady- state 
probleas,  in  which  the  cur\'e  C,.  is  merely  the  translation  of  with  uniform 

U  O 

velocity  along  the  .\-axis.  Then  S(z,t)  =  V^t  +  S(z-Vpt,0)  ([2, (8. 8)]),  so 
S.(z,t)  =  V  [1-S  (z-V  t,0)]  and  the  integration  in  (3.5)  and  (3.6)  can  be 
performed  e.\plicitly.  Such  problems  can  be  reduced  to  examination  of  S(z;0); 
they  can  also  be  formulated  directly  in  a  more  elementary  way  ( [13] ) . 

Examples  include  the  motion  and  shape  of  bubbles  and  fingers  in  llele-Shaw 
cells  of  finite  width. 

Problems  of  two- phase  flow  can  be  formulated  in  the  same  manner,  although 
the  Cauchy  data  on  the  interface  are  not  given  explicitly  in  terms  of  S,  since 
there  is  coupling  between  the  solutions  in  the  two  fluids. 
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Moving  boundary  problems  arise  also  in  many  other  contexts;  see,  for 
example,  the  introduction  to  [6] .  Some  of  these  involve  synthesis  or 
control,  and  the  present  approach  should  be  useful  there. 

Direct  problems  are  more  difficult  to  examine  in  the  present  formulation, 
for  they  require  the  determination  of  S  from  prescribed  sources-  Some 
consideration  has  been  given  to  this  general  problem  ( [6] ,  [9] ,  [10] ,  [15] , 
[16]).  Much  remains  to  be  done,  and  it  is  hoped  that  others  will  be 
encouraged  to  study  the  question . 
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A.B.  PARRY 

The  scattering  of  velocity  fields  by  an  airfoil 
in  compressible  flow 

1.  INTRODUCTION 

The  interaction  of  an  airfoil  with  a  fluctuating  velocity  field,  and  the 
calculation  of  the  resultant  scattered  field,  is  generally  of  great  importance 
in  both  aerodynamics  and  hydrodynamics.  Most  of  the  work  in  these  ares  has 
been  concentrated  on  low  frequency  interaction  problems  in  incompressible  or 
weakly  compressible  flow  (von  Karman  &  Sears  1938;  Sears  1940;  Kemp  1952, 
1973;  Osborne  1973;  Goldstein  1976,  chap.  3).  The  latter  can,  of  course, 
be  reduced  to  an  equivalent  incompressible  flow  problem  by  a  suitable 
transformation  (see,  for  example.  Ward  1955;  Landau  &  Lifshitz  1959).  Here, 
however,  we  consider  such  scattering  problems  in  as  much  as  they  apply  to 
noise  generation  by  blade  row  interactions  on  the  new  generation  of  advanced 
propellers.  In  this  application  the  blades  operate  in  the  high  subsonic 
regime  -  even  at  'take  off  and  'approach'  conditions.*  Accordingly,  the 
usual  methods  of  dealing  with  these  interaction  problems  are  inappropriate. 
Moreover  -  as  we  will  see  below  -  the  disturbance  velocity  field  is  not  always 
convected  with  the  mean  flow,  as  is  usually  the  case. 

The  way  in  which  the  fluctuating  velocity  field  is  modelled,  and  the 
far- field  sound  obtained  from  the  unsteady  pressure  distribution  across  the 
airfoil  chord  (and,  indeed,  across  the  whole  span  of  the  blade)  on  advanced 
propellers,  is  described  in  detail  in  Parry  (1988)  and  Parry  k  Crighton 

*  In  virtually  all  of  the  present  designs,  the  blade  tips  only  operate 
supersonically  at  the  'cruise'  condition  or  design  point. 
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(1989).  Here,  therefore,  it  is  only  necessary  to  consider  the  interaction  of 
a  single  Fourier  component  of  the  velocity  field  with  a  flat  plate.  The 
problem  is  assumed  to  be  two- dimensional. 2  Our  aim  is  to  calculate  the 
scattered  field  and  thereby  the  distribution  of  unsteady  pressure  across  the 
plate. 

In  §2  we  consider  the  interaction  of  an  airfoil  with  a  convected  gust  from 
upstream  infinity.  Applications  of  this  to  the  wake  interactions  on  a 
counter-rotation  propeller  show  considerable  differences  between  the  measured 
and  predicted  acoustic  field.  In  §3  the  work  is  extended  to  include  the 
interaction  of  an  airfoil  with  potential  (i.e.  non- convected)  velocity  fields 
from  both  upstream  and  downstream  infinity.  In  addition,  an  iterative 
technique  is  presented  which  can  be  used  to  put  the  solution  in  the  form  of  an 
asymptotic  series.  Comparisons  between  predictions  and  measurements  of  the 
far- field  noise  of  a  counter-rotation  propeller  show  that  the  predictions  are 
extremely  accurate,  in  terms  of  both  absolute  levels  and  noise  directivity. 

2.  CONVECTED  GUST  INTERACTIONS 

The  major  source  of  aerodynamic  interference  between  blade  rows  is  taken, 
usually  (and  naturally  enough),  to  be  the  unsteady  velocity  field  associated 
with  the  viscous  wakes  generated  by  the  upstream  blade  row. 3  We  suppose  that 


2  Some  justification  for  the  two-dimensional  approximation  comes  from  the 
asymptotic  analysis  of  Parry  k  Crighton  (1989)  who  showed  that  noise 
generation  is  highly  localised  at  discrete  radii. 

3  In  addition  the  tip- vortex,  of  current  interest  with  regard  to  advanced 
propellers,  can  be  represented  by  convected  gusts.  Tip  vortex  interactions 
can,  therefore,  also  be  described  by  the  methods  of  this  section  -  providing, 
of  course,  that  the  velocity  field  is  known. 
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the  gusts  -  representing  .the  velocity  field  -  are  of  relatively  high 
frequency;  this  is  appropriate  for  advanced  propellers  which  have  large 
numbers  of  blades  (in  addition,  we  must  remember  that  the  ear  is  most 
sensitive  to  the  higher  frequency  interactions  which,  therefore,  assume 
greater  importance). 

A  transformation 

We  start  with  a  single  Fourier  component  of  the  convected  velocity  field  given 
by 

V  =  u  exp[ik(Ut-x)]  (1) 

on  y  =  0.  The  axes  x  and  y  are  oriented  as  shown  in  figure  1  and  centred  on 
the  airfoil  leading  edge  with  U  the  velocity  in  the  x  direction  and  k  the 
wavenumber.  The  disturbance  velocity  potential  ^  will  satisfy  the  convected 
wave  equation  which  we  write  as 

9  9 

(1-M^)  -9  +  -0  -  2ikM^  -  ^  ,  2„2,  .  .  (2) 

df  ^  +  k  M  -  0  V  ; 

where  M  =  U/c^  is  the  Mach  number  and  the  time  dependence  e^^^^  is  implied. 

The  boundary  condition  on  the  airfoil  surface  is  simply 

^+v  =  0ony  =  0,  0<x<c  (3) 

where  c  is  the  chord  length  of  the  airfoil. 
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upstream  downstream 

Figure  1.  The  coordinate  axes  for  downstream  (leading  edge)  interactions  or 
upstream  (trailing  edge)  interactions. 


Ve  now  introduce  the  Prandtl- Glauert  coordinates,  scaled  on  the  semi- chord 
c/2,  given  by 

. _  X  =  2x/c,  Y  =  2i?y/c,  (4) 

where  ji  =  v(l-M^)  is  the  usual  compressibility  factor.  Ve  also  introduce  a 
new  potential  il>  such  that 


CU  r  am  ■) 

^(x,y)  =  2^  |>(X,Y)exp  i  — ^  ^  •  (5) 

t  1-  M 

The  result  of  (4)  and  (5)  is  that  (2)  then  reduces  to  the  standard  Helmholtz 
equation 

=  0  (6) 

where  the  new  nondimens ional  wavenumber  K  is  given  by 


cr 
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(7) 


K  = 


and  a  =  kc/2  is  the  reduced  frequency.  In  the  new  coordinate  system  the 
boundary  condition  on  the  plate  is  given  by 

=  -exp(-i/cX)  on  Y  =  0,  0  <  X  <  2,  (8) 

and  the  pressure  p  =  -pD^/Dt  =  -pU(ik+3/ftc)^  is  related  to  the  transformed 
potential  ii  by 


where  the  wavenumber  k  in  (8)  and  (9)  is 

a 


K  = 


(10) 


Airfoil  response 

The  solution  to  (6),  subject  to  the  boundary  condition  (8),  at  high 
frequencies  produces  a  pressure  distribution  which  oscillates  rapidly  away 
from  the  leading  edge  of  the  airfoil  where,  indeed,  it  will  be  (integrably) 
singular.  The  pressure  will  be,  therefore,  to  a  large  part  self- cancelling. 
The  trailing- edge  region,  where  a  Kutta  condition  is  applied,  should  then  be 
relatively  unimportant  so  that  the  pressure  distribution  is  much  the  same  as 
that  on  a  flat  plate  extending  to  downstream  infinity.  This  leading- edge 
problem,  or  two-part  boundary  problem,  can  be  solved  by  the  Wiener- Ilopf 
technique  (see,  for  example.  Noble  1958;  Grighton  1977).  The  solution  has, 
however,  been  given  previously  by  Landahl  (1961)  and  Goldstein  (1976),  who 
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used  an  alternative  approach,  and  we  simply  quote  the  result  as** 

2/>Uu  r  X  <tM 

Ap(X)  = 


[t(t(1+M)X]^/^ 


exp 


-^4 


^  1+M  ^ 


Heasurenent  vs.  prediction 

It  is  difficult  to  validate  this  result  directly  as  we  have  no  suitable  data 
on  blade  unsteady  response.  However,  by  combining  it  with  the  noise 
radiation  formulae  for  counter-rotation  propellers  (Hanson  1985),  and  an 
appropriate  model  for  the  (viscous  wake)  unsteady  velocity  field  (Parry  1988), 
we  can  obtain  predictions  of  far- field  noise.  (The  way  in  which  these 
different  stages  can  be  combined  into  a  robust  prediction  scheme  has  been 
described  by  Parry  k  Crighton  1989).  These  predictions  will  be  compared  with 
measurements  taken  from  the  flyover  tests  carried  out  by  Rolls-Royce  on  the 
Fairey  Gannet  counter-rotation  propeller:  details  of  these  tests  have  been 
discussed  by  Bradley  (1986).  This  comparison  has  been  given  before  by  Parry 
k  Crighton  (1989)  but,  for  the  sake  of  completeness,  and  since  the  comparison 
serves  as  a  check  on  the  blade  response  calculation,  we  will  give  it  again 
here  along  with  a  brief  discussion.  Since  the  front  and  rear  4- blade  rows  on 
the  Gannet  were  run  at  slightly  different  speeds,  the  interaction  tone 
components  could  all  be  separated  out  in  terms  of  frequency,  thus  allowing  us 
to  examine  each  tone  individually. 

The  first  interaction  tone  generated  by  the  Gannet  is  the  (1,1) 


Our  solution  is  the  complex  conjugate  of  Goldstein's  since  he  chose  a  time 

dependence  and  we  have  used  In  addition,  a  phase  term  e^^  is 

missing  from  our  result  since  here  the  reference  is  the  airfoil  leading  edge 
and  not  the  mid- chord  as  in  Goldstein's  work. 
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interaction. 5  The  directivity  plot  (sound  pressure  level  versus  radiation 
angle  where  0°  represents  the  direction  of  flight  or  the  propeller  axis)  for 
this  tone  is  shown  in  figure  2.  We  can  see  that  there  is  a  null  at  the  90® 
radiation  angle  in  both  the  measured  and  predicted  data.  This  is  to  be 

expected  since  the  (1,1)  interaction  tone  generates  a  plane- wave  mode  (on  the 
Gannet  which  has  equal  blade  numbers,  Bj  and  B2  on  the  front  and  rear  rows), 
i.e.  n^Bj  -  n2B2  =  0.  This  tone  peaks  on  the  propeller  axis  and  is  zero  in 
the  plane  of  the  rotor.  However,  predictions  are  typically  25dB  below  the 
measurements. 


Figure  2.  Gannet  measurements  vs  predicted  wake  interaction  noise  for  the 
(1,1)  interaction  tone. 


5  We  will  use  the  notation  (n^n^)  to  indicate  an  aerodynamic  interaction 
tone  generated  at  the  combination  frequency  njf,  +  n2f2  where  f j  and  f2  are 
the  blade  passing  frequencies  of  the  front  and  rear  blade  rows  respectively. 
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{  The  next  two  interaction  tones  generated  by  the  Gannet  are  the  (2,1)  and 

\ 

‘  (1,2)  interactions,  for  which  directivity  plots  are  shown  in  figures  3  and  4 

i 

'  respectively.  For  these  interaction  tones  the  predictions  are  at  least  lOdB 

below  the  measurements.  In  fact,  it  is  not  only  the  levels  which  are 
I  incorrectly  predicted  but  also  the  directivities:  predictions  are  20dB  below 

the  measurements  for  the  (1,2)  interaction  in  the  forward  arc,  and  40dB  below 
i  the  measurements  for  the  (2,1)  interaction  in  the  rear  arc. 


(2,1)  interaction  tone. 


134 


1 


Vc  BOV  refer  to  t&c  ^iscessloai  ob  iaieractiens  t<aie  directivities  'ivea 
%  Bradl<^'  (1SS6)  aad  Parry  k  Crigjktoo  (19§9).  l^ere  it  vas  t^t, 
altliOG^  tiic  neasared  directivities  arc  sisilar  ia  level  ia  forward  aad  rear 
arcs,  any  aoisc  soiree  oa  tfec  rear  blade  rov  Kill  prodacc  oaly  asynmetrical 
directivities  (except  for  tke  plaae  wave  case  vAcrc  Mj  =  02).  Tfccre  mast, 
tlicrcfors.  be  a  sigpificaat  aoisc  soiree  oa  tie  forward  blade  ror.  fa  order 
to  generate  toecs  at  tbc  (apa.^)  ioteractioa  toac  freqacacies,  tbc  noise 
scarce  can  only  be  that  dec  to  tbc  interaction  of  tbe  forvard  blade  roK  vitb 
the  potential  field  generated  by  tbc  rear  row- 

3.  POTEbTIAL  FIELD  IKrEIACriO.VS 

Accordingly,  vc  proceed  to  consider  the  scattering  of  a  potential  velocity 
field  by  an  airfoil.  In  addition  to  tbc  interaction  of  the  forward  blade  rev 
with  the  potential  field  generated  by  the  rear  row,  it  sccass  appiopriatc  to 
include  the  interaction  of  the  downstrean  blades  with  tbc  bound  potential 
field  generated  by  the  forward  blade  row,  because  that  is  likely  to  generate 
at  least  as  large  a  field  as  that  froa  the  interaction  of  the  forward  row  with 
the  potential  field  of  the  rear  row.  The  latter  involves  a  trailing  edge 
(weakly  loaded  if  a  Kutta  condition  is  satisfied),  the  fomer  a  leading  edge 
(highly  loaded).  As  before  we  assuac  that  the  detailed  velocity  field  is 
known  (the  modelling  of  the  bound  potential  flow  field  of  a  counter-rotation 
propeller,  in  c-onpressible  flow,  is  described  in  full  by  Parrj'  1988)  and 
consider  a  single  harmonic  component  of  the  upwash 

V  =  u  exp(iwt-i7x)  (12) 

on  y  =  0.  The  axes  x  and  y  are  centred  on  the  airfoil  trailing  edge,  for 
upstream  interactions,  or  on  the  leading  edge,  for  downstream  interactions,  as 
shown  in  figure  1.  The  velocity  in  the  x- direction  is  U  and  7  is  the  complex 


135 


vavcaas&cr. 


(1,2)  interaction  tone. 


Airfoil  response  -  opstrcaa  row 

An  analysis  similar  to  tliat  given  in  §2  shosws  that  the  transformed  velocity 
potential  related  to  ^  by  (5),  again  satisfies  the  Helmholtz  equation 

V  +  K  =  0.  The  boundary  condition,  however,  is  now  given  by 

* 

^  = -exp(-iK  X)  on  Y  =  0,  -ro  <  X  0  (13) 

where 

* 


■f 


[i 


(14) 


H  =  7c/2. 


For  sieplicity,  since  the  pressure  p  is  related  to  the  transformed  potential  ^ 
by  (9),  we  define  a  transformed  pressure 


The  solution  is  obtained  by  the  Wiener- Hopf  technique  and  will  be 
described  here  in  detail.  We  define  +  and  -  Fourier  transforms  by 

♦i(s,V)  =  r  |)(X,¥)H{a)e'«%  (16) 

where  H(X)  is  the  Heaviside  unit  function.  The  inverse  transform  is  given  by 

WS.V)  =  4  f  ♦(s,V)e>=*  ds.  (17) 

*'“00 


On  Fourier  transforming  the  Helmholtz  equation  (6),  we  obtain 

r'(s,v)  *  /»(s,Y)  =  0, 

where  we  use  primes  to  denote  differentiation  with  respect  to  Y  and 


(18) 

(19) 


Here  we  choose  the  branch  of  the  square  root  so  that  -ijsl  as  s  ±oo. 

The  branch  cuts  in  the  complex  plane  are  shown  in  figure  5.  The  wavenumber  K 
is  taken  to  have  a  small  imaginary  part. 
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Figure  5.  The  integration  contour  and  branch  cuts  in  the  complex  plane. 

We  now  consider  the  region  X  >  0  where  we  assume  continuity  of  pressure 

* 

across  the  wake.  From  (9)  and  (15)  p  is  also  continuous  across  the  wake  so 

that,  on  Fourier  transforming  (15),  we  obtain 

i(K-s)[f^(s,0-)  -  t^.(s,0+)]  =  0.  (20) 

Since  ii  must  be  an  odd  functin  of  Y  this  leads  to 

t^(s,0+j  =  0.  (21) 

* 

If  we  now  take  the  Fourier  transform  of  p  in  the  region  X  <  0  we  obtain 

-2i(K-s)f.(s,0+)  =  AP*(s)  (22) 

* 

where  we  have  again  used  the  fact  that  is  odd  in  Y  and  AP_  represents  the 

* 

transform  of  the  jump  in  p  across  the  airfoil. 

The  boundary  condition  (13)  can  be  Fourier  transformed  to  give 
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fO 

(s,0)  =- 

•'-00 


e 


S-/C 


(23) 


where  the  integral  converges  provided  that  .Sfs)  <  S(k  )  and  S(z)  denotes  the 
imaginary  part  of  z. 

From  (6)  and  the  requirement  that  the  fields  decay  at  infinity  we  obtain 


#(s,Y)  =  C(s)e 


rUY 


(24) 


where  we  have  taken  Y  to  be  positive.  (For  Y  negative  we  can  again  use  the 
fact  that  i>  is  odd  in  Y.)  On  differentiating  (24)  and  setting  Y  =  0+  we 
obtain 

♦  '(s,0)  =-i:t*(s,0+).  (25) 

Substituting  (21)  and  (23)  into  (25)  we  find  that 

♦;(s,0)  +-^=  i;t».(s,0+).  (26) 

S-K 

We  now  put  x  =  \/(K-fs) (K- s)  so  that,  on  dividing  through  by  /(K-s) ,  (26) 
becomes 

♦;(s,o)  i 


=  iy(K+s)f_ (s,0+) . 


(27) 


(s-/)V5^ 

Here  the  first  term  on  the  left-hand  side  is  a  +  function  and  the  right-hand 
side  is  a  -  function.  The  second  term  on  the  left-hand  side  can  be  split,  in 
the  usual  way,  into  the  sum  of  a  +  function  and  a  -  function.  We  thus 
rewrite  (27)  as 

*1  /"o  {\\ 

i  r  1  1 


♦;(s,0) 


y(M  (s-«  )  /(k./) 


=  -  iy(K+s)*_  (s,0+)  - 


(s-K 


(28) 
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where  the  left-hand  side  is  a  +  function  and  the  right-hand  side  is  a  - 
function.  By  the  usual  arguments  resulting  from  Liouville's  theorem,  both 
sides  of  (28)  are  equal  to  an  entire  function  E(s)  which  must  be  identically 
equal  to  zero  (otherwise  ♦|(s,0)  would  diverge  as  s  -»  oo  implying  that  ii  has  a 
singularity  at  the  origin).  From  (22)  and  the  right-hand  side  of  (28)  we 
then  obtain 


AP_(s)  = 


2i(/£-s) 


(s-/c  )(K+s) 

We  now  apply  the  inverse  Fourier  transform,  as  defined  by  (17),  whence 


Ap*(X) 


r. 


(k-s) 


isX 


/(k-/) 


ds. 


(29) 


(30) 


On  wrapping  the  integration  contour  around  the  branch  cut  in  the  upper  half 
plane  we  find  that 


Ap  (X)  = 


2giKX+iT/4  (K+K-is') 


'00 

Jo 


Xgi  _  2[k-k  ) 


ds'  - 


* 

(31) 


kJ  (K-/C  )  ^  v^(s'+iK+i/c  ) 
where  the  final  term  represents  the  contribution  from  the  pole  at  s  =  «;  .  On 

evaluating  the  integral  we  obtain 


Ap*(X)  =  ^  l  (e^^^wfi^^X|(iK^•k' 


J^K 


*2 


2giKX- ij/d 


vf 


(32) 


)r(K+K  )  |X|] 


where  w(x)  =  e'^  erfc  (-ix)  is  the  complex  error  functin  (see  Abramowitz  k 
Stegun  1965). 

We  now  impose  a  Kutta  condition  at  the  trailing  edge  of  the  airfoil.  In 
the  appendix  we  show  that  this  is  equivalent,  in  this  case,  to  removing  the 
inverse  square- root  singularity  at  the  trailing  edge.  Then,  from  (9),  the 
pressure  jump  across  the  airfoil  is  given  by 
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Ap(X)  = 


(33) 


where  we  have  defined  the  normalised  frequencies 

(TM  ~ 

'  °  {ny  ’  " '  it+mj  • 


(34) 


Leading  edge  correction 

In  the  previous  section  we  dealt  with  just  the  trailing  edge  problem  and 
consequently  the  leading  edge  effects  have  been  neglected.  Since  leading 
edge  effects  can  be  important,  as  they  are  usually  heavily  loaded,  we  will 
discuss  a  correction  to  the  present  results  to  account  for  such  effects. 

We  will  use  a  technique  developed  by  Landahl  (1961)  and  Adamczyk  (1974) , 
and  discussed  by  Amiet  (1975),  for  downstream  convected  gust  interactions. 
This  involves  an  iterative  technique  for  the  solution  of  a  3  part  boundary 
problem.  The  situation  is  shown  in  figure  6  with  the  upwash  specified  on 
-2  _<  X  _<  0. 
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Figure  6.  Three- part  boundary  value  problem. 


The  first  iteration  involves  solving  the  trailing  edge  problem  with  upwash 

specified  on  -<»  <  X  <  0;  this  is  the  case  discussed  in  the  previous  section, 

* 

i.e.  we  set  Ap^  =  Ap,5^/5Y  =  -exp(-i/c  X),  on  Y  =  0.  The  second  iteration 
involves  correcting  the  upstream  boundary  condition  on  -w  <  X  <  -2  without 
affecting  the  boundary  condition  on  the  airfoil,  i.e.  on-2<X<0.  We 
therefore  require  the  'new'  pressure  difference  across  -oo  <  X  <  -2  to  be  minus 
that  obtained  on  the  first  iteration,  i.e.  Ap2  =  -Ap^  on-oo<X<-2,  and  the 
upwash  on  -2  <  X  <  00  to  be  zero.^  This  second  iteration  produces  an  error  in 


6  We  will  use  the  suffix  2  throughout  this  section  to  denote  'second 
iteration'  values. 
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the  boundary  condition  on  0  <  X  <  «  which  could  be  corrected  by  a  third 

” 

i  iteration  and  so  on. 

f 

We  define  a  new  coordinate  X  centred  on  the  airfoil  leading  edge  so  that 


j 

i 


X  =  X  +  2.  The  pressure  jump  Ap2  across  -oo  <  X  <  0  is  minus  that  obtained  in 
the  previous  section  so  that 


Ap2(X)  =  -Ap(X)  =  -Ap(X-2).  (35) 

Then,  from  (9)  and  (15),  the  normalised  pressures  are  related  by  a  phase 

shift:  Ap2(X)  =  -exp[-2i(rM^/(l-M^)]Ap  (X-2).  From  (32)  this  leads  to^ 


Ap.(X) 


2{k-k  ) 

9  *9 

(K^-/c  ^) 


-2i 


at 


1-M" 


X  |e^'^(^"2)w[iv(x-2|(iK+i«*)]-e'^'^  (36) 

Since  we  are  considering  a  high  frequency  problem,  a  and  n  are  both  large. 
Then,  from  (7)  and  (14),  the  argument  of  the  complex  error  function  in  (36)  is 

also  large  on  X  <  0  so  that,  from  Abroraowitz  Ic  Stegun  (1965),  we  can  use  the 
approximation 

»[iix-2|(iRH-i/)]  ~  ■■■■■ . . . :  .  (37) 

\4iX-2| (iK+i«  ) 

(In  order  to  use  this  approximation  we  have  used  the  fact  that 

|arg  v(iK+iA;  )  |  <  x/2.)  In  addition,  since  {i  is  large  and  S(ii)  0,  the  last 

term  in  braces  in  (36)  is  from  (14),  exponentially  small  on  X  <  0; 

* 

consequently,  this  term  will  be  neglected.  The  jump  in  P2  across  X  <  0  can 


7  Recall  that  the  inverse  square  root  singularity  term  has  been  removed  in 
order  to  satisfy  the  Kutta  condition. 
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therefore  be  approximated  by 


4p!(X) 


2{k-k  ) 


exp 


/  ^  1“  M 

yT(iK+i/c  )(r-/c  ^)|X-2| 

The  boundary  condition  is  that  of  no  upwash  so  that 


■  ar  -  n 
-2i — ^+iK(X-2).  (38) 

1-M^  J 


—  =  0, 


X  >  0, 


The  Fourier  transform  Ap2(X)  is  given  by 


4P2-  (®)  =  - 


2{k-k  ) 


J 


r(iK+i«*)(K^-/^) 


exp 


-2i 


(tM" 


1-M" 


2iK 


rO  e 


i(s+K)X 


(39) 


where,  for  consistency,  we  have  replaced  Y  with  Y. 

The  problem  is  therefore  defined  by  the  Helmholtz  equation  with  the 

* 

boundary  condition  (39)  and  the  jump  in  P2  across  X  <  0  given  by  (38). 
solution  is  again  obtained  by  the  Wiener- Hopf  technique  with  +  and  - 
transforms  defined  by  (16)  and  the  inverse  transform  defined  by  (17). 


The 


-00 


dX. 


(40) 

We  now  use  the  fact  that  K  is  large  so  that  the  integrand  in  (40)  oscillates 
rapidly.  Then  the  integral  is  dominated  by  contributions  from  close  to  X  =  0 
and,  by  following  Murray  (1974),  we  can  approximate  (40)  to  leading  order  in 
l/(K+s),  by 


AP‘  (s) 


2{k-k  ) 


/ 


ir(iK+i/)(K^-/^)(s+K) 


aT 

exp[-2i[^  +  k]].  (41) 


From  this  point  we  proceed  with  the  Wiener- llopf  technique  in  the  usual 
manner  and  the  pressure  jump  Ap2(X)  is  finally  obtained  as 
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AP2(X)  =  - ^  ,  :::....:r=:r 

(i/i+i(r  )v  5r(l-M^)  (i<r- i/i) 


y 

* 

where  we  have  reintroduced  the  normalised  frequencies  a,  fi,  c  and  <r  and,  in 
addition,  we  have  defined 


a  = 


1-M' 


a  = 


(ItM) 


(43) 


We  note,  from  (42),  that  we  have  restored  the  inverse  square  root  singularity 
at  the  airfoil  leading  edge.  The  pressure  distribution  on  the  airfoil  is 
given  by  the  sum  of  (33)  and  (42).  Comparing  these  equations  we  see  that, 
apart  from  the  exponential  decay  term  e'^^^  in  p  and  the  inverse  square  root 
singularity  l/^^l+X)  in  P2,  the  first  and  second  stage  solutions  are  basically 
of  the  same  form  apart  from  the  factor  lfJ{in-i-i<T  )  in  P2.  Since  we  are 
considering  a  high  frequency  problem  both  //  and  a  (and  hence  a  )  are  large. 

In  the  high  frequency  limit  therefore  the  first  correction  to  the  trailing 
edge  problem  is  0(1/V^)  smaller  than  the  leading  order  term  -  even  though  the 
correction  term  includes  the  inverse  square  root  singularity  at  the  airfoil 
leading  edge.  We  conclude  that,  in  the  high  frequency  limit,  the 
semi- intinite  airfoil  model  is  a  valid  approximation  and  provides  accurate 
results  to  leading  order  in  a. 


Airfoil  response  -  downstream  row 

The  airfoil  response  calculation  for  the  downstream  blades  is  calculated  in 
much  the  same  way  as  for  the  upstream  blades.  The  only  difference  is  that 
now  we  have  to  solve  a  leading- edge  problem  instead  of  a  trailing- edge 
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problem.  The  least  singular  solution  is  chosen,  and  has  an  inverse  square 
root  pressure  singularity  at  the  leading  edge. 

The  Wiener- Hopf  technique  can  be  used,  as  before,  to  obtain  the  pressure 
difference  across  the  airfoil  (which  now  occupies  y  =  0,  0<x<oo)  as 


Measurement  vs.  prediction 

As  in  §2  we  combine  the  results  given  here  for  airfoil  response  with  the 
appropriate  noise  radiation  formulae  and  a  model  for  the  unsteady  potential 
velocity  field  (Hanson  1985;  Parry  1988)  in  order  to  obtain  predictions  of 
far-field  noise  which  can  be  compared  with  the  Gannet  data.  Note  that  here, 
unlike  the  case  considered  in  §2,  we  have  noise  sources  on  both  front  and 
rear  blade  rows.  The  relative  phasing  of  the  two  sources  should,  therefore, 
be  corrected  to  account  for  spatial  separation  (as  discussed  by  Hanson  1985). 
However  we  will,  for  the  present,  consider  the  two  fields  separately. 

The  first  interaction  tone  generated  by  the  Gannet  is  the  (1,1) 
interaction.  The  far-field  directivity  of  this  tone  is  shown  in  figure  7. 


interaction  noise  for  the  (1,1)  interaction  tone. 

The  first  thing  to  note  is  that  the  predicted  potential  field  interaction 
noise  is  significantly  greater  than  the  predicted  wake  interaction  noise:  the 
predicted  forward  and  rearward  potential  field  interaction  tones  are, 
typically,  20dB  greater  than  the  predicted  wake  interaction  tone.  The 
predicted  potential  field  interaction  noise  agrees  extremely  well  with  the 
measured  data  in  both  forward  and  rear  arcs  -  except,  perhaps,  for  a 
discrepancy  in  the  range  140  -  160  degrees. 

The  next  two  interaction  tones  generated  by  the  Gannet  are  the  (2,1)  and 
(1,2)  interaction  tones,  for  which  directivity  plots  are  shown  in  figures  8 
and  9.  Here  we  see,  again,  that  the  predicted  potential  field  interaction 
noise  levels  are  significantly  greater  than  the  predicted  wake  interaction 
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noise  levels  -  typically  by  20dB  (and  more  at  some  angles).  Note  that,  as  we 
remarked  in  §2,  sources  on  the  forward  blade  row  (due  to  the  upstream 
potential  field  interaction)  generate  far- field  directivities  different  from 
those  of  sources  on  the  rear  blade  row  (due  to  the  downstream  row  and 
potential  field  interactions). 

We  emphasise  that  the  predictions,  shown  in  figures  7-9,  are  absolute 
level  predictions  governed  solely  by  the  theoretical  prediction  scheme 
outlined  above  and  dependent  on  the  high-frequency  approximation.  The  inputs 
for  the  calculation  of  the  scattered  fields  are  the  incoming  velocity  fields: 
models  for  these  are  described  by  Parry  (1988). 
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interaction  noise  for  the  (2,1)  interaction  tone. 


ANGLE  0 
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Figure  9.  Gannet  measurcraents  vs.  predicted  wake  and  potential  field 
interaction  noise  for  the  (1,2)  interaction  tone. 

4.  CONCLUSIONS 

We  have  described  a  model  for  the  calculation  of  unsteady  velocity  fields 
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scattered  from  airfoils  in  subsonic  coapressible  flow  at  high  reduced 
frequencies.  The  aodel  covers  interactions  of  leading  edges  or  trailing 
edges  with  convected  or  nonconvected  (potential)  gusts.  In  addition,  an 
iterative  technique  has  been  outlined  by  which  the  scattered  field  can  be  put 
in  the  fora  of  an  asyaptotic  series  with  successive  terns  decreasing  by 
O(lVff),  where  <r  is  the  reduced  frequency. 

Coaparison  with  noise  aeasureaents  taken  from  a  Fairey- Gannet 
counter-rotation  propeller  has  shown  that  the  analysis  produces  extremely 
accurate  results  -  in  terms  of  both  the  absolute  level  and  the  far- field 
directivities  -  with  no  adjustment  whatsoever  of  the  theoretical  predictions. 

In  addition,  the  comparison  with  measured  data  has  shown  that,  for  the 
Gannet,  the  downstream  wake  does  not  dominate  the  aerodynamic  interactions, 
and  that  the  potential  flow  field  around  each  row  generates  significant  and 
indeed  dominant  effects,  both  upstream  and  downstream. 
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APPENDIX:  THE  KUTTA  CONDITION 

In  our  analysis  the  airfoils  have  been  modelled  as  flat  plates  (finite  or 

semi- inf inite) .  As  we  have  found,  the  pressure  jump  across  the  plate  has 

inverse  square  root  singularities  at  the  edges  of  the  plate.  In  order  to 

alleviate  the  singularity  at  the  trailing  edge  we  must  introduce  a  vortex 

sheet,  extending  to  downstream  infinity,  across  which  the  tangential 

velocities  jump  but  the  pressure  is  continuous.  If  the  strength  of  the 

vortex  sheet  is  fixed  in  order  to  cancel  exactly  the  trailing  edge 

singularity,  then  a  Kutta  condition  is  said  to  be  satisfied. 8  The  use  of  a 

Kutta  condition  in  unsteady  flow  is  a  matter  of  controversy  at  the  present. 

For  the  moment,  however,  we  will  assume  that  a  Kutta  condition  is  satisfied. 

We  look  for  a  potential  which  satisfies  the  Helmholtz  equation  (6)and 

is  odd  in  Y.  Then,  downstream  of  the  trailing  edge,  there  is  a  jump  in  the 

potential  across  the  vortex  wake  so  that 

ii/»X 

=  ±Ge  ^  on  Y  =  0±,  X  >  0,  (Al) 

* 

where  and  G  are  to  be  determined.  Now  p^(X)  is  continuous  across  Y  =  0 
since  Pj^(X)  is  continuous^  so  that  (15)  implies  that  =  -k.  Since  there  is 

8  A  more  detailed  discussion  of  the  Kutta  condition  in  unsteady  flow  is 
provided  by  Crighton  (1981,  1985). 

* 

8  The  definitions  of  pj^  and  Pj^  are  the  same  as  in  §§2,  3  except  that  we  have 

introduced  a  subscript  K  on  those  parameters  relating  to  the  velocity 
potential 
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no  additional  upwash  the  boundary  condition  is 


•^y—  =  0  on  Y  =  0,  X<0.  (A2) 

We  now  have,  once  again,  a  two-part  boundary  value  problem  which  we  will 
solve  by  using  the  Wiener- Hopf  technique. We  then  obtain,  on  X  <  0, 

This  shows  that  the  effect  of  the  Kutta  condition  here  is  merely  to  remove 
the  inverse  square  root  singularity  at  the  trailing  edge,  i.e.  we  select 

G  =  -  (A4) 

V(K+/c)(K-«  ) 

so  that  the  sum  of  (32)  and  (A3)  contains  no  term  in  l/Tfxi. 


A.  B.  Parry 

Department  of  Mathematics 
University  of  Strathclyde 
Glasgow  G1  IXII 
Scotland,  U.K. 


10  Crighton  (1977,  Chap.  91  shows  how  the  Wiener- llopf  technique  can  be  used 
to  solve  a  trailing  edge  proolem  in  unsteady  compressible  flow  with  a  Kutta 
condition  imposed. 
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M.K.  PIDCOCK 

Boundary  problems  in  electrical  impedance 
tomography 

ABSTRACT 

The  reconstruction  problem  for  Electrical  Impedance  Tomography  is  an  extremely 
ill- posed  non-linear  inverse  problem  and  the  results  obtained  are  highly 
sensitive  to  modelling  and  measurement  errors  associated  with  the  technique. 

In  this  paper  we  continue  our  investigations  into  some  simple  problems  which 
involve  geometric  errors  and  which  can  be  solved  by  using  perturbation  theory. 
The  sensitivity  of  the  technique  to  these  errors  is  made  explicit  and  we 
describe  one  way  of  overcoming  these  difficulties  which  is  suggested  by  the 
analysis. 

1.  INTRODUCTION 

Electrical  Impedance  Tomography  (EIT)  is  a  technique  of  medical  imaging  which 
uses  the  contrast  in  the  electrical  conductivity  of  different  body  tissues  to 
produce  an  image  of  the  conductivity  distribution  within  a  part  of  the  body. 

In  many  circumstances  this  image  can  be  interpreted  as  a  physical  image  and  a 
number  of  potential  medical  applications  of  this  technique  are  being 
investigated.  The  data  used  to  obtain  these  images  are  measurements  taken  on 
the  surface  of  an  object  of  the  electrical  potential  which  are  induced  in  the 
object  by  the  application  of  known  electrical  currents  to  that  surface. 

Mathematically,  the  problem  of  BIT  can  be  posed  in  the  following  way. 
Suppose  that  an  object  fi  with  boundary  Xl  consists  of  an  isotropic  Ohmic 
material  with  conductivity  distribution  a.  If  (5  is  the  electrical  potential 
in  U  then 
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where 


V.(ffV^)  =  0 


in  it 

^  is  known  (measured)  on  Xl 

-(T  ^  =  j  is  known  (applied)  on  dO. 

It  is  important  to  note  that  one  consequence  of  this  model  is  that  the 
applied  current,  j,  satisfies 

j  ds  =  0.  (1) 

The  aim  of  EIT  is  to  determine  a  from  the  boundary  data.  It  has  been 
shown  [1-5]  that  this  problem  has  a  solution  for  many  reasonable  conductivity 
distributions  and  reconstruction  algorithms  have  been  described  [6-11]  which 
work  well  in  simple  cases  and  in  the  absence  of  errors.  However,  the 
reconstruction  problem  of  EIT  is  a  highly  non-linear  inverse  problem  and  it 
has  been  demonstrated  numerically  [7]  that  the  Fr^chet  derivative  of  the 
mapping  which  takes  an  applied  current  to  a  measured  voltage  has  singular 
values  which  decay  exponentially.  The  inversion  procedure  is  therefore 
extremely  sensitive  to  the  possible  errors  associated  with  this  technique. 

There  are  a  wide  range  of  such  errors  which  have  to  be  considered.  These 
range  from  basic  deficiences  in  the  mathematical  model  of  the  system  to 
numerical  errors  introduced  in  the  solution  of  equations  in  the  reconstruction 
algorithm  and  data  measurement  errors  defined  by  the  instrument  specification. 
Each  has  its  own  characteristics  and  a  detailed  study  of  them  is  essential  if 
EIT  is  to  become  a  useful  diagnostic  tool.  The  interaction  of  any 
reconstruction  technique  with  a  model  of  these  errors  represents  one  of  the 
attractions  of  inverse  problems. 

In  this  paper  we  will  consider  just  one  type  of  error  -  that  of  imperfect 
knowledge  of  boundary  shape.  We  will  consider  a  class  of  simple  problems 


156 


which  can  be  solved  using  perturbation  methods.  This  may  be  a  little 
artificial  but  it  does  enable  us  to  follow  a  parameter  which  characterises  the 
boundary  error  in  an  explicit  way.  In  EIT  the  electric  current  is  applied  to 
the  surface  via  a  number  of  electrodes  positioned  on  the  boundary  but  in  this 
work  we  will  assume  that  the  current  is  adjustable  at  any  point  on  the 
boundary. 

In  Section  2  we  will  describe  the  basic  problem  and  the  ideas  behind  our 
calculations  and  in  Section  3  we  will  discuss  a  particular  example  of  boundary 
error  where  the  boundary  is  perturbed  from  its  assumed  shape  in  a  simple  way. 
In  Section  4  we  will  extend  our  analysis  to  more  general  perturbations  and 
suggest  a  scheme  for  identifying  the  parameters  in  the  perturbation.  We  use 
polar  coordinates  (r,0)  throughout. 


2.  THE  MODEL  PROBLEM 

Most  of  the  studies  in  EIT  have  been  concerned  with  the  two-dimensional 
problem  and  a  simple  example  often  considered  is  that  of  distinguishing 
between  two  conductivity  distributions  disc,  0  <  r  <  1 

defined  by 

(rj(r,0)  =  1  0  <  r  <  1 

(r2(r,^)  =  ff  0<r<R<l 

1  R  <  r  <  1 

It  is  easy  to  show  that  if  the  applied  current  j(0)  is  given  by 


00 


j(0)  =  S  a  cos(n^)  +  b  sin(n0) 
n=l  “ 


then 


00  r"f 


— |^ajjCOs(n^)  +  bj^sin(n<?) 
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and 


«S(r,^;ir2,j)  =  Ajjr^[aj^cos(n0)  +  bjjSin(n^)] 


where 


1 

\  =  s 


ri-/,RV"i 

with  s  =  r  (r>R), 

and  11 

R  (r<R) 

<r-l 

(T+l 


As  we  are  able  to  choose  coefficients  (a^^,  b^^)  in  anyway,  it  is  natural  to 
ask  if  there  are  any  combinations  which  are  somehow  better  than  others. 

Gisser  et  al  [12]  have  suggested  that  if  we  wish  to  distinguish  between  and 
then  the  best  normalised  currents  (i.e.  llj||  =  1)  to  use  are  those  which 
maximise 


where 


1  p 

<f>g>  =  T  fg  ds  and 


iifi|2  =  <  i,f>. 


In  other  words,  we  should  try  to  maximise  some  average  difference  between 
the  signals  measured  in  the  two  cases.  Other  criteria  for  the  choice  of 
optimal  currents  have  been  proposed  [13]  and  these  will  have  different 
stability  properties  when  interacting  with  the  geometric  errors  considered  in 
this  paper.  We  have  not  yet  studied  these  alternative  currents  and  we 
confine  our  attention  to  those  proposed  in  [12] . 

For  the  case  of  the  problem  described  above,  these  optimal  currents  turn 
out  to  be  trigonometric  functions  =  cos(m^),  sin(m0)  for  integer  m,  and 

2 

the  corresponding  value  of  6  =6-  is  - - iy—.  It  is  interesting  to  note 

^  (l+/iR^"’) 

that  if  this  value  of  6^  is  less  than  the  accuracy,  E,  of  the  measuring 
equipment,  then  and  a 2  •'^t  distinguishable  using  the  current  since 

the  value  of  could  be  an  effect  of  noise.  However,  if  >  E  then  we  can 
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distinguish  between  these  two  distributions.  This  is  the  basis  of  the 
calculations  which  we  will  report  in  the  next  section.  We  will  calculate  6 

m 

arising  when  these  optimal  trigonometric  currents  are  applied  to  a  body  whose 
shape  is  thought  to  be  the  unit  disc  and  which  consists  of  material  with  a 
uniform  conductivity  distribution,  If  the  boundary  were  correctly  known 

1 

and  all  the  measurements  were  accurate  then  the  value  of  .6  would  be  zero. 

m 

s  The  error  in  the  boundary  description  will  mean  that  a  non- zero  value  of 
,  will  be  obtained  and  if  5  >  E  we  could,  erroneously,  infer  that  there  is  a 

f  ^ 

;  circular  anomaly  of  the  form 

I 

3.  BOUNDARY  DISTORTION 

Consider  the  situation  where  the  angular  displacement  of  the  drive  electrodes 
is  correct  but  the  polar  description  of  the  boundary  curve  is  r  =  t{9)  rather 
'  than  r  =  1.  In  this  case  we  have  that 

1  f  t'(^) 

^  ^  (2) 

and 

■  [  jds  =  w(0)j(^)  d^ 

^  '*0 

f  0  O 

where  v{$)  =  r{6)  +  r'^(0) 

If  r{6)  =  1  +  eF(^)  then  we  can  write  (2)  in  the  form 

d<^  2  S 

^(l+eF(^),^)  =  +  ftje  +  026  +  +  ... 

where  0^,0^...  are  functions  of  6  and  the  partial  derivatives  of  with 

respect  to  r  and  6,  evaluated  at  the  point  (l+eF(0),^). 

We  can  estimate  the  partial  derivatives  at  this  point  by  expanding  these 

functions  in  a  Taylor  series  about  the  point  (1,^).  If  we  then  write  ^{r,6) 
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in  the  form 

^{r,0)  =  f{T,e)  +  /(r,^)e  +  +  ... 

we  find,  finally,  that 

^(l+eF(^),^)  =  Pq  +  Pi^  +  ^2^^  ^3^^  ■*■... 

where  the  functions  P^,P^,...  are  complicated  functions  of  6  and  the  partial 
derivatives  of  with  respect  to  r  and  9  evaluated  at  the  point  (1,^). 

Suppose  that  we  have  the  simple  case  t{6)  =  1  +  e  cos(0),  i.e. 

F(0)  =  cos(0),  then  we  find  that 
f 

j  ds  =  1+2€  cos  d+e^  j(0)  d^. 

J^n  -^0  *■ 

If  we  now  try  to  apply  the  optimal  currents  described  earlier,  we 
find  that  the  condition  (1)  can  be  satisfied  exactly  only  for  =  sin(m0), 
m  an  integer.  So,  for  example,  if  we  apply  jj(^)  =  sin(0)  to  the  distorted 
boundary  we  find  that  at  the  point  (1,9) 

df 

^0  =  ^ 

5/  d<f> 

^1  =  ^  +  cos  ^  ^  +  sin(^)  ^ 
or 

and  similar,  but  more  complicated,  expressions  for  other  P's. 

Since  ^  satisfies  Laplace's  equation  in  fl  so  too  do  .  On  the 

boundary,  r  =  1,  we  find  that  P^  =  sm{9)  and  P^  ~  P2  ~  After 

considerable  computation  it  follows  that 

=  ||^^(l+e  cos(9),9)  -  cos(0)|| 

e  89 

=  4(1+  9  €  +  ...) 

In  general,  if  j(0)  =  sin(m0)  then  we  find  that 
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l^m  “ 


1 

f— 1 

Im^-ll 

[  2j 

1/2 


c  +  0(e'^) 


1 


We  can  generalise  the  type  of  perturbation  slightly  by  considering 
boundaries  defined  by  r(^)  =  1  +  e  cos(k^)  (k  =  1,2,...).  If  we  try  to  apply 
the  optimal  currents  we  find  that  equation  (1)  restricts  these  currents 

to  sin(m^)  for  m  an  integer  and  that 


k  =  4  * 


m=k 


,  2,2  1/2 

k  rin  +k  ^  o 


in>k 


1  C  1  2^1/2  q 

priyj  ]  e  +  0(e  ) 


(m+k) 


m<k 


In  Table  1  we  give  values  of  for  a  range  of  values  of  k  and  m.  It  is 
clear  that  there  is  a  direct  relationship  between  the  detectable  error  in 
boundary  shape  and  the  measurement  accuracy  of  the  system.  It  is  interesting 
to  note  that  the  observed  behaviour  gives  a  possible  scheme  for  identifying 
displacements  of  this  rather  special  type.  It  appears  that  if  we  apply 
currents  of  increasing  spatial  frequency  all  we  need  to  do  is  to  locate  the 
dip  in  the  values  of  A  against  frequency  in  order  to  identify  k.  Further 
comments  on  the  interpretation  of  j^A^  can  be  found  in  Pidcock  and  Breckon 
[14]. 
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Q| 

1 

2 

3 

4 

5 

6 

1 

0.25 

0.53 

0.28 

0.19 

0.15 

0.12 

2 

3.02 

0.25 

1.02 

0.53 

0.36 

0.28 

3 

3.77 

4.51 

0.25 

1.52 

0.77 

0.53 

4 

4.68 

5.01 

6.01 

0.25 

2.01 

1.02 

5 

5.64 

5.84 

6.26 

7.51 

0.25 

2.51 

6 

6.61 

6.76 

7.01 

7.51 

9.00 

0.25 

Table  1.  Values  of  k  for  various  values  of  k  anu  ;ii.  All  entries  should 


be  multiplied  by  e. 

4.  GENERAL  FOURIER  PERTURBATIONS 

It  is  interesting  to  note  that  for  small  e,  the  relationship  between  and  e 
given  the  previous  section  is  essentially  linear  and  that  it  appears  possible 
to  perform  simple  experiments  to  determine  the  parameters  of  the  perturbation 
once  its  general  form  is  known.  The  relative  ease  of  these  tests  encourages 
us  to  think  that  more  general  perturbations  which  can  be  described  in  terms  of 
a  Fourier  expansion  might  be  identified  by  a  suitable  series  of  such  tests. 
Such  an  expansion  should  be  very  appropriate  in  the  case  of  che  human  body 
where  the  relative  smoothness  of  the  body  surface  should  lead  to  an  economical 
description  in  terms  of  trigonometric  series. 

Consider,  therefore  the  more  general  boundary  perturbation  given  by 

1  00 

^  A  cos(n0) 
n=l  " 

If  we  try  to  apply  one  of  the  optimal  boundary  currents  it  is  easy  to  see 
that  (1)  implies,  once  again  that  we  must  use  sin(m^)  for  integer  m. 

Following  the  analysis  described  in  the  previous  section  to  first  order  in  e 
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ve  fled  tiat 


^9 


8^  8^ 

h  =  ^  wCi’ 

at 


aad  it  folloKS  tS^t 


9{fS  =  5-sw(«^ 


=  0 


d^ir.0)  =  -  E  ^  Bfsja(W). 
k=l  ^  ^ 

vbcre  4((‘-‘>  Vtr  <‘"'5  M- 

Hcncc,  the  pctcatial  actually  Bcasurcd  differs  friws  that  expected  *f  there 
vere  no  boundary  error  by  an  aaonnt 

Sojid)  =  ^{U(S{B),e)  ^  e  g^{ut?{e),e)  -  ^  o(£^) 


=  c  E 
k-1 


c°  sin(k^)  -f  O(c^) 


(3) 


«here  cj  =  -  {2k^l)  A,,,.]. 

Ve  now  have  a  possible  schene  to  determine  the  Fourier  coefficients.  It 
goes  as  follows-  Apply  a  series  of  currents  sin(m^)  to  the  object  and 
measure  resulting  voltages  at  a  number  of  points-  Use  a  Discrete  Fourier 
Transform  to  express  the  measured  voltage  in  terms  of  its  Fourier 

components  and  use  the  above  expression  (3)  as  a  system  of  linear  equations  to 
determine  {Aj^}.  The  implementation  of  this  scheme  and  a  detailed  study  of 
its  numerical  stability,  together  with  an  investigation  into  the  effects  of 
using  only  a  finite  number  of  electrodes, is  the  subject  of  future  work  in  this 
area. 
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OHCUiSlBE 

Ve  iiftve  scca  that  the  rccoastrectioa  problca  for  BIT  is  cxtrcsely  seasitivc  to 

errors  ia  t&e  boaadary  slaape.  Ilic  aaalysls  prcseatcd  bas,  kHrevcr.  offered  a 

possible  vay  to  overcooe  this  problea.  1  stady  of  tbis  netbod  vill  be 

preseated  elsevbere. 
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B.D.  SLEEMAN 

Interior  and  exterior  inverse  problems  for 
the  Helmholtz  equation 

1.  namucriM 

In  this  paper  we  discuss  a  nuaber  of  recent  developacnts  relating  to  inverse 
prouleas  for  the  Helaholtz  equation.  In  particular  we  concentrate  on  the 
problcB  of  detentining  the  geonetry  of  an  unknown  doaain  (e.g.  vibrating 
Bcabrane  or  scattering  obstacle)  fro«  given  data. 

The  paper  is  presented  in  two  parts.  In  Part  I  we  consider  the  classic 
inverse  problca  of  determining  an  unknown  doaain  fi  froa  a  knowledge  of  the 
eigenvalues  of  the  laplacian  defined  in  Q.  In  §1  we  survey  the  classic 
asysptotic  estimates  for  the  counting  function  N(I)-  Beginning  with  tlie 
fundamental  results  of  Hermann  Vcyl  we  survey  the  most  recent  results  for  non 
smooth  domains.  In  §2  we  take  up  the  Veyl- Berry  conjecture  regarding  the 
asymptotics  of  N(A)  for  fractal  domains.  In  particular  we  describe  the 
recent  contributions  of  Fleckinger  and  Lapidus. 

Part  II  of  the  paper  is  concerned  with  the  important  inverse  acoustic 
scattering  problem.  In  §3  we  formulate  the  direct  scattering  problem  which 
provides  the  basic  setting  for  the  inverse  problem.  §4  discusses  the  central 
question  of  uniqueness  of  reconstruction  of  an  unknown  scattering  obstacle 
from  far  field  data. 

In  §5  we  consider  the  "exterior"  analogue  of  the  asymptotics  of  N(A)  by 
discussing  the  high  frequency  behaviour  of  the  scattering  phase  s(k).  In 
particular  we  concentrate  on  the  asymptotics  of  s(k)  for  non- smooth  domains, 
which  complements  the  recent  results  of  Melrose,  and  also  fractal  domains. 

This  latter  result  provides  the  exterior  analogue  of  the  Veyl- Berry 
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conjecture.  §6  returns  to  the  practical  problca  of  devising  algorithas  for 
the  nuecrical  reconstruction  of  the  unknown  scattering  obstacle. 

Vhile  this  paper  is  largely  e.xpository,  a  nuaber  of  the  results  are  new 
and  have  not  previously  appeared  in  the  literature.  Theorem  2.3  is  new  as  is 
Theorem  4.4  regarding  uniqueness.  The  results  embodied  in  Theorems  (5.2), 
(5.3)  and  (5.4)  are  also  uew. 


PART  I 

INTERIOR  INVERSE  PROBLEMS 
§1  Geometry  of  the  Counting  Function 

Let  Q  be  an  arbitrary  non  empty  bounded  open  connected  set  in  R"(n  >  1)  with 
boundary  F  =  59  and  consider  the  eigenvalue  problem 

-Au  =  Au  in  9, 

u  =  0  on  r,  (1.1) 

wlierc  As  T,  d  /dxt  denotes  the  Dirichlet  Laplacian  in  9.  The  parameter  A  is 
k=l 

said  to  be  an  eigenvalue  of  the  problem  (1.1)  if  there  exists  a  u  0  in 
satisfying  -Au  =  Au  in  the  distribution  sense.  It  is  well  known  that  the 
spectrum  of  (1.1)  is  discrete  and  consists  of  an  infinite  sequence  of 
eigenvalues  which  may  be  ordered  according  to  their  multiplicity  as 


0  <  Aj  <  A2  -  <  <  • •• 


where  A^  h  «  as  i  -»  00. 


In  1912  Veyl  [47,48]  estaoiished  the  classical  result  that 

2/n 


Aj~C„ 


where  =  (2t)^(Bjj)  depends  only  on  the  dimension  'n'.  Here  |9|j^ 
denotes  the  n- dimensional  Lebesgue  measure  or  "volume"  of  9  and  is  the 


w 


as  1  00, 


(1.2) 
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voluac  of  the  unit  ball  in  R°. 

Another  way  of  estiaating  the  asyaptotic  behaviour  of  the  eigenvalues 
is  to  work  with  the  “counting  function"  N(A)  defined  as 

NO)  =  #{i  >  1.0  <  for  A  >  0.  (1.3) 

Consequently  Veyl's  result  can  be  restated  as 

N(A)  ~  (2x)‘”  as  A  -  00.  (1.4) 

In  a  classic  paper  entitled  "Can  one  hear  the  shape  of  a  drum?"  Kac  [19] 
has  asked  the  following  question:  Can  soaeone  with  perfect  pitch  recover  the 
precise  shape  of  a  drum  just  by  listening  to  its  fundamental  tone  and  all  the 
overtones?  This  question  has  motivated  some  important  advances  in  the  last 
two  decades.  In  the  first  place,  it  is  natural  to  ask  whether  the  problem 
has  a  unique  solution.  Unfortunately,  the  answer  appears  to  be  no  in 
general.  Urakawa  [44]  has  discovered  two  isospectral  domains  in  R"  (n  >  4) 
which  are  not  isometric.  Despite  this  it  is  possible  to  recover  a  lot  of 
topological  information  about  Q  from  the  spectrum  of  (l.T)  and  in  particular 
from  the  counting  function  and  other  related  functions.  Indeed,  if  F  is 
smooth  (i.e.  of  class  C”)  then  Seeley  [36]  and  Pham  The  Lai  [24]  have  shown 
that 

N(A)  =  (2x)'"  +  0(a("'^)/^)  as  A  h  a>.  (1.5) 

The  proof  of  this  result  makes  use  of  techniques  from  the  theory  of  spectral 
transforms  and  of  Fourier  integral  operators  (c.f.  Hormander  [16]).  More 
recently  Ivrii  [17]  has  shown  that  if  fl  is  a  bounded  domain  with  C“  boundary  F 

and  if  ft  does  not  have  too  many  multiply  reflected  closed  geodesics  then 

N(A)  =  (2t)'"  BJftJA"/2-C;|F|jj.jA("'^)/2+o(A("‘^)/2  as  A  -»  oo  (1.6) 
where  is  a  positive  constant  depending  only  on  'n'.  There  are  a  number  of 
extensions  of  this  result.  For  example  Ivrii 's  result  extends  to  the  Neumann 
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problem,  the  only  difference  being  that  the  minus  sign  of  the  second  term  on 
the  right  hand  side  of  (1.6)  is  changed  to  a  plus  sign.  Results  are  also 
known  for  the  "impedance"  or  Robin  boundary  value  problem  [41]  as  well  as  for 
higher  order  positive  elliptic  operators  with  locally  constant  leading 
coefficients  [12]. 

An  alternative  attack  on  the  problem  is  to  study  the  asymptotics  as  t  0” 
of  the  "partition  function"  (or  trace  of  the  heat  semigroup) 


z(t)  =  r 
Jo 


■It  ®  -A-t 

dN(A)  =  S  e  ^ 
i=l 


(1.7) 


which,  when  the  integral  e-xists,  may  provide  more  information  than  N(A). 

Thus  for  example  if  F  is  C“  and  n  =  2,  McKean  and  Singer  [31]  and  others  have 
shown  that 

1  00 


i®l  1^1  A  00  .  /« 

*  sd-")  *  j!,  Ci‘  ^ 


as  t  -»  0 


(1.8) 


here  h  is  the  connectivity  of  fl  and  the  coefficients  are  metric  invariants. 
Indeed  they  are  polynomials  in  the  curvature  of  F  and  its  derivatives. 

Smith  [42]  has  shown  how  to  compute  these  coefficients  using  symbolic 
manipulation  techniques. 

If  F  is  not  smooth  then  neither  (1.6)  or  (1.8)  are  expected  to  hold.  For 
example  if  fl  has  an  outward  pointing  cusp  then  Vaechter  [46]  has  shown  that 
(1.8)  takes  the  form 

|ft|  |r| 


m  =  477  - 


+  0(t'''),  0  <  u  <  1/2. 


(l.Q) 


8(Tt)V2 

If  we  consider  the  Neumann  problem  then  even  the  first  term  of  (1.5)  may  not 

hold  if  the  boundary  F  is  "too  long".  This  can  be  demonstrated  in  the 

following  example  due  to  Fleckinger  and  Metivier  [11] . 

For  a  given  positive  number  /?  define  the  set 

^/?  =  {(x,y)  e  K^'|x  6  (0,l),0<y<l  +  2  j'\  (x)|,  (1.10) 

'  ^  j6N  f 
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where  (Ij)j  G  N  is  an  infinite  sequence  of  disjoint  open  intervals  in  (0,1) 
and  where  denotes  the  characteristic  function  for  the  set  I-,  i.e. 

j  ^ 

rl,  X  6  I. 

'’ij  °  (O,  X  (!  Ij 

then  for  the  Neumann  problem  ve  find  that  for  0  <fi  <  1/2 

N(A)  s;  as  A  -» 00  (1-11) 

where  a  means  that  there  exist  two  positive  constants  G'  and  C"  such  that 

<  N(A)  < 

for  all  A  sufficiently  large. 

In  the  following  section  we  discuss  the  asymptotics  of  N(A)  and  Z(t)  in 
the  extreme  case  when  F  is  "fractal”. 

2.  TDE  WEYI^ BERRY  CONJECTURE  AND  FRACTAL  DOMAINS 

In  1979  M.  V.  Berry  [3]  motivated  by  the  study  of  the  scattering  of  light 
by  random  surfaces  conjectured  that  if  F  is  "fractal"  with  Hausdorff  dimension 
II  G  (n-l,n)  then  (1.6)  takes  the  form 

N(A)  =  (2T)'"BJft|jjA"/2  .  Gjj^„//j,(F)a”/2  +  o(a"/2)  as  A  h  oo,  (2.1) 
here  is  a  positive  constant,  depending  only  on  n  and  H,  and  //jj(F)  denotes 
the  II- dimensional  Hausdorff  measure  of  F. 

Note  that  if  F  is  sufficiently  smooth  then  II  =  n  -  1  and  (2.1)  reduces  to 
(1.6).  However  in  ge.neral  F  is  very  irregular  and  hence  11  >  n  -1.  Thus 
Berry's  conjecture  seems  a  reasonable  one.  However  in  (1986)  Brossard  and 
Carmona  [4],  through  a  series  of  illuminating  examples,  showed  that  (2.i) 
cannot  be  true  in  general  and  proposed  that  the  Hausdorff  measure  (dimension) 
should  be  replaced  by  the  less  familiar  Minkowski  measure  (dimension  S) . 

In  order  to  understand  recent  contributions  to  the  Weyl- Berry  conjecture 
we  define  the  Hausdorff  and  Minkowski  dimension  as  follows 
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Definition  2.1 


{00 

inf  E  r®L  where  the  infimum  is  taken  over 

£-40' 

all  countable  coverings  of  A  by  open  balls  of  radius  r-  <  e.  The 

number  /tjj(A)  in  [0,oo)  is  called  the  d- dimensional  Hausdorff  measure  of  A  and 
1!(A)  =  inf{d  >  0,/f^(A)  =  0}  =  sup  {d  >  0,^j(A)  =  oo}  is  called  the  Hausdorff 
Dimension  of  A. 

Definition  2.2 

Given  e  <  0  let  =  {x  G  R”:d(x,r)  <  t}  be  the  e- neighbourhood  of  T. 

For  d  >  0  let 

=  li«  (2.2) 

€-40'^ 

be  the  d- dimensional  upper  Minkowski  content  of  T.  Then 

D  =  D(r)  =  inf{d  >  0,/tj(r)  =  0} 

=  sup{d  >0,/fj(r)  =  oo} 
is  called  the  Minkowski  dimension  of  F. 

If  0  <  /fjj  <  00  then  r  is  said  to  be  Minkowski  measurable  and  is  called 
the  Minkowski  measure  of  F. 

Examples 

1.  Let  A  be  the  set  of  rational  numbers  in  [0.1],  then  H(A)  =  0  and  D(A)  =  1. 


2.  Let  A  be  the  set  A  =  n  K.,  where  K-  is  the  union  of  2^  disjoint  intervals 

i=0  ^  ^ 

of  length  ttj  such  that  =  1  and  <  a^/2.  Then 

llfA'l  =  lim  inf  r-TTT-Y-  =  lim 

If  we  have  the  classic  1/3- Cantor  set  then  =  3  ^  and 
11(A)  =  D(A)  =  log  2/log  3. 


3-  Lapidus  [25] 

Let  a  be  an  arbitrary  fixed  positive  number  and  let  ft  C  R"  be  the  bounded 
open  set 


and 


Then 


a  =  U  (I,  «  J)  where  J  =  n  >  1 

1=1  * 


li 


=  {{W) 


-a 


=  n  -  1  and  D(ft)  =  n  -  1  +  (a+l)"^. 


4.  Brossard- Carmona  [4] 

2 

Let  ft  c  R  be  the  countable  disjoint  union  of  all  the  small  open  cubes 

belonging  to  the  successive  generations  defined  as  follows 

Let  be  a  nondecreasing  sequence  of  positive  integers.  The  0-th 

generation  contains  1  square  of  side  1.  The  1-st  generation  contains  4 

squares  each  of  side  1/3  and  is  divided  into  (Pj)  congruent  small  squares. 

Similarly  the  i-th  generation  consists  of  4x5^’^  squares  of  side  3'^  and  is 
2 

divided  into  congruent  small  squares  and  so  on. 

Brossard  and  Carmona  [4]  show  that  irrespective  of  the  sequence  {P^^} 

II (ft)  =  log  5/log  3.  However  if  ,  for  example,  P^  =  [a^]  for  same  a  >  1  and 

f 

a 

where  [  ]  indicates  "integer  part  of"  then  D(ft)  =  log  5a  /log  3a. 

It  is  known  and  of  course  clear  from  the  above  examples  that  11(A)  <  D(A). 
More  than  this,  examples  3  and  4  show  that  while  11(A)  is  constant  D(A)  is 
parameter  dependent.  The  significance  of  this  is  that  the  Minkowski 
dimension  is  more  sensitive  to  the  "roughness"  of  the  boundary  of  ft  than  is 
the  llausdorff  dimension. 

Returning  now  to  our  main  theme  regarding  the  Weyl- Berry  conjecture 
lapidus  and  Fleckinger-Pelle  [26]  have  proved 
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TheoreB  2.1 

If  ft  is  bounded  and  if  dit  is  fractal  with  Minkowski  dimension  6  G  (n-l,n) 

then 

N(A)  =  (2r)"''  Bj^lftljjA"/^  +  0(A^/^)  as  A  oo. 

In  addition  this  result  has  been  extended  by  Lapidus  [25]  to  more  general 
elliptic  operators  and  to  the  Neumann  boundary  value  problem.  More  recently 
Fleckinger  [10]  has  proved. 

Theorem  2.2 

If  ft  is  bounded  and  if  5ft  is  5-Minkowski  measurable  with  5- Minkowski 
measure  /i  then 

IN(A)  -  (2t)'"  Bj^|ft|j^A"/^|  <  C(n,5)/iA^/^  for  all  A  >  A^. 

In  order  to  give  the  reader  an  idea  of  the  arguments  used  to  establish 

results  such  as  Theorems  2.1  and  2.2  we  outline  the  proof  of  the  folloying 

result  concerning  the  asymptotics  of  the  partition  function  Z(t)  for  fractal 
2 

domains  in  [R  .  This  result  is  analagous  to  Theorem  2.2. 

Theorem  2.3 

2 

If  ft  C  R  is  bounded  and  its  boundary  T  is  5- Minkowski  measurable  with 
5- Minkowski  measure  n,  then  there  exists  a  constant  7  depending  on  6  so  that 

1^1  A /9 

z(t)  =  ^  <  7/it"®'  for  all  t  <  t^. 

In  order  to  prove  the  theorem  we  need  the  following  preliminary  results. 
(1)  Dirichlet- Neumann  Bracketing 

2 

Suppose  ft  is  a  bounded  open  set  in  R  ,  let  NQ(A,-A,ft)  be  the  counting 
function  for  the  Dirichlet  Laplacian  on  ft  and  Nj(A,-A,ft)  th  counting  function 
for  the  Neumann  Laplacian  on  ft.  Then  it  is  well  known  (c.f .  Gourant  and 
Hilbert  [9])  that  the  following  proposition  holds. 


Proposition  2.1 


If  ft'  c  ft  then 

N^(A,-A,ft)  >  N^(A,-A,ft').  (2.3) 


If  ftj  and  ft2  are  two  disjoint  open  sets  in  ft  with  ft  =  ft^  U  ft2  then 
NQ(A,-A,ftj)  +  NQ(A,-A,ft2)  <  N^(A,-A,ft) 

<  Nj(A,-A,ft)  (2.4) 

<  Nj(A,-A,ftj)  +  Nj(A,-A,ft2) 

(2)  Polygonal  Domains 
Proposition  2.2  [45] 

Let 


(i)  D  be  a  polygonal  domain  with  boundary  5D 

(ii)  Pp...,Pj^  be  the  vertices  of  and  let  be  the  infinite  wedge  of 

angle  7^  with  vertex  P^  such  that  the  boundary  of  the  wedge  contains 
the  two  edges  adjacent  to  P^. 

(iii)  Define  for  y  >  0,  7  =  min  7^ 

Bj(y)  =  {A  e  V.|(l(A,Pj)  <  y) 

1  n 

R  =  -  sup{y|B.(y)  n  B.(y)  =  for  all  i  it  j,  U  B.  (y)  c  D) 

2  J  k=l 

Then 


2  2 
n  ^-'^i 


(R  sin  7/2)' 


|D|  \dd\  n  "  ’^i  1D|  1  «  -  (0 

■  4irt  ^  I4ir^  |<(5n+20  x  jgt  • 

Note  [5]  that  a  similar  but  less  precise  result  holds  for  the  "Neumann" 
partition  function  together  with  a  change  of  sign  in  the  boundary  term  of 
(2.6). 


Outline  of  the  Proof  of  Theorem  2.3 

To  begin  with  we  introduce  a  positive  number  such  that  for  all 
e  e  (0,6^) 
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<  2ii 


(2.6) 


-P. 


aiid  choose  Pq  G  N  such  that  2  ^  In  (2.6)  we  define 

=  {x  G  R^|d(x,«2)  <  e} 

where  d(.,^)  denotes  the  Euclidean  distance  to  the  boundary  5ft. 

2 

For  each  integer  p  we  consider  a  tessellation  of  R  into  congruent 
non- overlapping  squares 

2  .  .  - (P"Po) 


(2.7) 


Q>  jC  ^  2  side  =  2 
Sp  P  P 


Define 


*0  = 


%  =■ 


Aj  =  {(i .  C  =  %  u  [  u 


ftj  =  ft\ftj, 


ICpSAp  '■pj 


ft"  =  ft\ft^ 
P  '  P 


We  also  define  the  boundary  sets 

Bp  =  {Cp  e  Z^IQ/  n  ea  t  o,g^  n  =  a), 
R.  =  u  Q,  . 

"  (p€Bp<p 

We  now  make  the  following  observations  and  estimates 


01)  ft"  C  ft  with  €  =  ^/^r]  =  ^.2 
r  P  P 

Furthermore  R„  c  ft, 

P 


-  (P+Po) 


where  ft  is  defined  by  (2.7) 
P 


175 


